UNCLASSIFIED 

SECURITY  CLASSIFICATION  Of  THIS  PAGE 


/  r: 

~"/V- 


14.  REPORT  SECURITY  C*  OSSIFICATION 

UNCLASSIFIED 


2a.  SECURITY  CLASSIFICATION  AUTHORITY 


REPORT  DOCUMENTATION  PAGE 


lb.  RESTRICTIVE  MARKINGS 


3.  DISTRIBUTION' AVAILABILITY  OF  REPORT 


2b.  DECLASSIFICATION  DOWNGRADING  SCHEDULE 


Approved  for  public  release; 
distribution  unlimited. 


4.  PERFORMING  ORGANIZATION  REPORT  NUMBER(S) 

NSWC  TR  87-27 


6a.  NAME  OF  PERFORMING  ORGANIZATION  6b.  OFr!CE  SYMBOL 

(If  applicable) 

Naval  Surface  Weapons  Center  K104 


6c  ADDRESS  (City.  State ,  and  ZIP  Code) 

Dahlgren,  Virginia  22448-5000 


5.  MONITORING  ORGANIZATION  REPORT  NUMBER(S) 


7a.  NAME  OF  MONITORING  ORGANIZATION 


7b.  ADDRESS  (City.  State,  and  ZIP  Code) 


8 a  NAME  OF  FUNDING  SPONSORING  gb.  OFFICE  SYMBOL 

ORGANIZATION  (If  applicable) 

Naval  Surface  Weapons  Center  K104 


9.  PROCUREMENT  INSTRUMENT  IDENTIFICATION  NUMBER 


8c  ADDRESS  (City,  State,  and  ZIP  Code) 

Dahlgren,  Virginia  22448-5000 


10.  SOURCE  OF  FUNDING  NOS 


PROGRAM 
ELEMENT  NO. 


PROJECT 

TASK 

NO 

NO- 

WORK  UNIT 
NO. 


1  TITLE  (Include  Security  Classification) 

INTEGRATION  OF  THE  TRIVARIATE  NORMAL  DISTRIBUTION  OVER  AN  OFFSET  SPHERE  AND  AN 

12.  PERSONAL  AUTHOR(S) 

Armido  R.  DiDonato 

13a.  TYPE  OF  REPORT 

13b.  TIME  COVERED 

14  OATE  OF  REPORT (Yr..  Mo..  Oay ) 

IS.  PAGE  COUNT 

Final 

FROM  TO 

1988,  February 

91 

16.  SUPPLEMENTARY  NOTATION 


COSATI  CODES 


18.  SUBJECT  TERMS  (Continue  on  reverse  if  necessary  and  identify  by  block  number) 


/"I  r  .... 


19.  A8STRACT  (Continue  on  reverse  if  necessary  and  identify  by  block  number)  ! 

;  A  numerical  procedure  is  given  for  evaluating  the  integral  of  the  trivariate 
normal  distribution  over  an  offset  sphere  of  radius  R.  A  numerical  algorithm  is 
also  described  for  computing  R  when  P  is  given.  The  procedures  generally  give  P 
or  R  to  six  decimal  digits.  Program  listings  in  HP  BASIC,  for  the  HP  9845B 
computer,  are  included.  -  ,  i  .  .  %  f  -i 


oT'Si 


■lecte 

NJG2  2®® 


20.  DISTRIBUTION  AVAILABILITY  OF  ABSTRACT 
UNCLASSIFIED  UNLIMITED  SAME  AS  RPT.  OTIC  USERS 
22a.  NAME  OF  RESPONSIBLE  INDIVIDUAL 


21.  ABSTRACT  SECURITY  CLASSIFICATION 

UNCLASSIFIED 


22b  TELEPHONE  NUMBER 
(Include  Area  Code ) 


I  A.  R.  DiDonato 

DO  FORM  1473, 84  MAR 


(703)  663-8036 

EDITION  OF  i  APR  83  15  OBSOLETE 


lie  OFFICE  SYMBOL 

K104 


UNCLASSIFIED 

SECURITY  CLASSIFICATION  OF  this  PAGE 


Best 

Available 

Copy 


.**V* , 


(.••j'av'Jljfcg' 


va.a,*:#Vj 


I  *>4  .'l  .'i  .V|  J|  .*>  * •****! k*  j|‘ , 


NSWC  TR  87-27 


FOREWORD 


The  work  described  in  this  report  was  done  in  the  Space  and  Surface  Systems 
Division  of  the  Strategic  Systems  Department.  Its  purpose  was  to  design  software, 
suitable  for  a  high  quality  mathematics  and/or  statistical  subroutine  library, 
which  yields  probabilities  and  percentage  points  for  3-dimensional  normal 
distributions  over  off-set  spheres. 

The  author  is  grateful  to  Mel  Holloway  for  his  guidance  and  direct  assistance 
with  the  technical  word  processor  that  was  used  to  produce  the  draft  copy  of  this 
report . 
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I.  INTRODUCTION 

In  the  application  of  the  methods  of  operation  research  to  problems  of 
military  strategy,  it  is  sometimes  necessary  to  compute  the  value  of  a  kill 
probability  expressible  in  the  form 


p(m)  F (x,y , z,u , v,w)  dx  dy  dz. 


where 


IF(x,y,z,u,v,w)  =  - 7 -  E(x/u)  E(y/v)  E(z/w) 

2tt  VJtT  uvw  (2) 

E(t)  =  exp(-t2/2). 

Here  P  is  the  kill  probability.  The  shots  or  bursts  are  assumed  to  form  an 
uncorrelated  trivariate  normal  distribution  with  their  mean  or  expected  burst 
point  set  at  the  origin  of  a  rectangular  cartesian  coordinate  system  Oxyz.  The 
distribution  has  standard  deviations  u,  v,  w  in  the  x,  y,  z  directions,  respec¬ 
tively.  This  configuration  may  be  assumed  without  loss  of  generality  for  any 
trivariate  normal  distribution,  since  if  the  distribution  is  correlated  and  does 
not  have  its  mean  at  the  origin,  a  translation  and  a  rotation  of  axes  will  bring 
about  the  situation  which  has  been  described. 

A  target  T,  which  may  be  either  a  point  or  area  target,  is  assumed  in  some 
arbitrary  position  in  the  Oxyz  space.  The  distance  m  is  the  radial  or  slant 
range  miss  distance  of  an  arbitrary  burst  point  (x,y,z)  from  the  target.  If 
a  target  is  a  point  at  a  position  (H,K,L),  then  m  is  the  ordinary  distance 
[(x  -  H)2+  (y  -  K)2+  (z  -  D2]1/2.  if  t  is  an  area  target,  m  is  the  minimum 
distance  from  (x,y,z)  to  points  of  the  target.  For  area  targets  m  is  set  to 
zero  for  all  burst  points  (x,y,z)  lying  within  the  target. 

The  expression  p(m)  is  a  function  of  m  as  prescribed  on  the  basis  of  physical 
conditions  or  assumptions  of  the  problem.  Possible  forms  of  this  function  which 
have  practical  applicability  are:  (1)  p(m)  =1  if  m  <  ft  (a  constant),  p(m)  =  0 
otherwise;  (2)  p(m)  is  a  decreasing  linear  or  exponential  function  of  m;  (3)  p(m) 
is  a  monotonically  decreasing  step  function. 
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This  report  describes  a  procedure  for  the  efficient  computation  of  these  kill 
probabilities  by  a  high  speed  digital  computer  for  an  important  practical  case, 
that  in  which  the  target  is  a  point  at  an  arbitrary  position  (H,K,L)  where, 
without  loss  of  generality,  H,  K,  L  are  nonnegative  and  the  function  p(m)  has  the 
first  form  mentioned  above,  i.e.,  p(m)  =  0  for  m  >  R,  p(m)  =  1  for  m  <  R.  The 
constant  R  may  be  thought  of  as  the  lethal  radius  of  the  weapon. 

The  function  p(m)  has  the  effect  of  reducing  the  field  of  integration  to  the 
interior  of  the  sphere  S  of  radius  R  with  the  target  (H,K,L)  as  center,  i.e., 

S  :  (x  -  H)2  +  (y  -  K)2  +  (z  -  L)2  =  R2  .  (3) 

The  kill  probability  P  is  then  given  by 
/•L+R  /-  k+Y  /*H+X 

P  =  |  I  L  F(x,y,z,u,v,w)  dx  dy  dz,  (4) 

'  L-R-J  K-Y  •'H-X 

where 


=-  vV  - 


(y  -  K)2-  (z  -  L)2, 


(5) 


Rather  than  carry  out  the  numerical  triple  integration  of  Equation  4 
directly  we  make  use  of  an  available  computer  program  that  yields  the  probability 
Pc(r,H,K,u,v)  of  a  shot  falling  under  a  bivariate  normal  distribution  inside  a 
circle  in  the  Oxy-plane  of  radius  r  and  center  (H,K).  Pc  is  the  two-dimensional 
analog  of  P.^»^»^. 

Geometrically,  one  observes  then  that  P  can  be  obtained  by  considering 
circular  slices  of  S  parallel  to  the  Oxy-plane.  For  a  fixed  z  in  [L  -  R,  L  +  R], 
the  xy-integration  in  Equation  4  over  a  slice  of  radius  r  =  \/r?  -  (z  -  L)2  yields 
Pc.  Weighting  Pc  with  1/(V  2tt  w)  E(z/w)  and  integrating  the  result  over  z  in 
[L  -  R,  L+R]  gives  P,  i.e., 

I  fL+R 

P=  -  l  E(z/w)Pc(r,H,K,u,v)  dz,  (6) 

\FT*  wJb-R 

where 


2 
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.  c(L+R)  /V~2  w 

-  I  ECVTt)  Pc(r,H,K,u,v)  dt, 

VrT  •'(L-R)/V,~2  w  r 

[From  (2),  E(V2t) 


=  exp(-t2)] 


r  =VR2  -  (L  -  V2wt)2. 


Hereafter,  we  assume  Pc  is  available  and  the  evaluation  of  P  reduces  to  the 


single  integration  indicated  in  Equation  6  or  Equation  8. 

The  symmetry  properties  of  F  (see  Equation  2)  allow  H,  K,  L  to  always  be 
taken  as  nonnegative.  Since  F  >  0  and  bounded,  the  order  of  integration  in 
Equation  h  is  immaterial.  Thus,  as  long  as  H  is  associated  with  u,  K  with  v, 
and  L  with  w,  it  does  not  matter  thereafter  which  is  called  L  and  w.  For  example, 
if  the  order  of  integration  is  chosen  so  that  the  x-integration  is  performed  last, 
then  if  initially  H  =  10,  u  =  5,  L  =  20,  w  =  7,  we  simply  let  H  =  20,  u  =  7, 

L  =  10,  w  =  5.  In  this  way,  we  can  always  refer  to  Equation  6  or  Equation  8  as 
the  basic  representation  for  P,  where  L  and  w  are  always  associated  with  the  z-  or 
t- integration  indicated  in  Equation  6  or  Equation  8  with  the  understanding  that 
the  order  of  integration  may  have  been  changed  and  the  variables  H,  K,  L  and  u,  v, 
w  renamed  as  in  the  example  above. 


II.  COMPUTATION  OF  P 


Sometimes  we  shall  write  P(R)  or  P(R,H,K,L,u,v,w)  for  P.  As  stated  in  the 
previous  section,  R  denotes  the  radius  of  the  sphere  S  given  by  Equation  3.  The 
trivariate  normal  distribution  F  is  defined  in  Equation  2  by  its  mean  (0,0,0)  and 
its  standard  deviations  u,v,w  in  the  x,  y,  and  z  directions,  respectively.  By 
appropriate  normalization  of  the  variables  R,  H,  K,  L,  u,  v,  w  it  is  easy  to  show 
that  P  is  in  general  a  function  of  six  independent  variables.  The  program  for 
computing  P  is  called  ELLC0V. 

In  most  cases  P  will  be  computed  by  the  numerical  integration  of  Equation  8, 
where  it  is  assumed  Pc  is  available.  There  are  two  special  cases  however  where 


P  can  be  evaluated  in  closed  form.  We  shall  call  these  cases  A  and  B.  The  deri¬ 


vations  of  the  results  for  these  will  be  given  in  Appendix  A.  Case  A  is  found  in 
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where 


daws  x  =  E 


rx 

(vTt) 

j  o 


exp(t2)  dt  (Dawson's  integral  [1;  p.  298],  [3]) 


jerfc  x  =  2/%/tT  f  exp(-t2 

'  x 

L2  =  v(l-  (w/u)2 | 

L4  =  (L  +  R)/(V?w) 

L5  =  (L  -  R)/(vTw) 

So  =  L/  (\/2  wL2  )  -  RL2/(n/2w) 
Fi  =  L/(V2wL2)  +  RL2/(V2w) 


exp(-t2)  dt  =  1  -  erf  x. 


(19.1) 

(19.2) 

(19.3) 

(19.4) 

(19.5) 


In  order  to  achieve  improved  accuracy,  Equations  10-19,  are  used  in  modified  form 
in  the  computing  program,  ELLCOV.  It  is  also  not  difficult  to  show  that  Equa¬ 
tions  13,  14,  16  approach  Equation  10  as  u  +  w. 

If  cases  A  and  B  do  not  apply  an  attempt  is  made  to  improve  efficiency  by 
sensing  if  Pc  is  numerically  constant  over  the  effective  range  of  the  t-integra- 
tion  in  Equation  8.  Call  this  case  C.  We  discuss  this  case  below  after  specify¬ 
ing  the  effective  range  of  the  t-integration ,  (see  Equation  29). 

In  general,  excluding  cases  A,  B  and  C,  P  is  computed  by  numerical  integra¬ 
tion.  The  24-point  Gaussian  quadrature  rule  Cl,  p.  916]  is  used  for  this  purpose 
with,  in  some  cases,  an  additional  24-point  Gaussian. 

The  most  time  consuming  part  of  computing  the  integrand  in  Equation  8  is 
the  evaluation  of  Pc,  the  elliptic  coverage  function^*^.  jn  most  cases  Pc  is 
obtained  from  the  subprogram  function  ELLCV.  In  the  special  situations  that 
u=vorH=K=0,  (referring  to  Equation  8)  Pc  is  called  the  circular  coverage 
function,  GIRGV.  and  the  generalized  circular  error  function,  GCEF,  respective- 
ly4,5,6.  Both  are  computed  by  the  subprogram  CIRCV.  In  these  special  situations, 
a  parameter  V  is  set  to  1  and  0,  respectively,  otherwise  V  =  2.  When  V  =  0  or  1 , 
Pc  can  be  computed  about  10  times  faster  than  the  case  when  V  =  2.  Hence,  since 
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the  order  of  integration  in  Equation  8  can  be  changed,  the  appropriate  inter¬ 
changes  are  always  made  to  make  V  =  1,  if  possible.  If  this  cannot  be  done 
(u  y  v,  v  y  w,  w  y  u)  then  an  attempt  is  made  to  obtain  V  =  0.  However  some 
instability  can  appear  in  this  case,  when  computing  GCEF,  if 

1 (-5 )  >  u/v  >  1(5).  (1 (N )  =  10N)  (20) 

Hence  if  V  =  0  and  Inequality  20  holds  then  the  computation  of  Pc  is  carried  out 
as  described  below  for  the  case  V  =  2. 

If  neither  V  =  1  nor  V  =  0  can  be  brought  about  then  the  order  of  integration 
is  determined  by  finding  the  largest  ratio  a/b  <  1,  where  a,  b  take  the  values  of 
u,  v,  w.  The  pair  from  u,  v,  w  giving  the  largest  ratio  (<1)  is  then  used  in  the 
first  two  integrations  (the  evaluation  of  Pc )  with  the  corresponding  H,  K,  L 
values.  For  example  if  u  =  2,  H  =  0,  v  =  5,  K  =  2,  w  =  3,  L  =  6,  then  the  a/b 
ratios  are  2/5,  3/5,  2/3.  Hence  v  and  w,  K  and  L,  are  interchanged,  the  integrand 
in  Equation  8  becomes 

exp(VTt)  P  (r.0,6,2,3) ,  (21) 

where 

r  =  ;  (22) 

the  limits  of  integration  become  (2  -  R)/5\/~2,  (2  +  K)/by/~2.  If  the  two  largest 
ratios  are  equal,  then  the  one  with  the  largest  denominator  is  used  to  fix  the 
first  two  integrations. 

Choosing  the  integration  order  in  this  way  permits  the  most  accurate  computa¬ 
tion  of  Pc.  This  remark  is  based  on  an  empirical  study  in  Reference  4,  where  it 
was  noted  that  the  numerical  quadrature  of  Pc  is  best  when  the  standard  deviations 
are  nearly  equal  and  not  small.  We  mention  at  this  point  that  the  program  ELLCV 
has  been  up-graded  over  its  description  in  Reference  4.  Namely,  the  accuracy  for 
Pc  has  been  increased  from  6  to  approximately  8-decimal  digits.  This  was  required 
in  order  to  deal  effectively  with  the  inverse  problem  discussed  in  Section  III. 

It  is  possible,  in  some  cases,  by  the  tests  given  below  to  determine  if 
P  <  2.5(-8)  or  P  >  1  -  2.5(-8).  In  such  situations  P  is  set  to  zero  and  one, 
respectively.  The  tests  below  are  discussed  in  Appendix  B. 

6 


P  set  to  0  if  R2  <  D2  and  R"  exp(-ct2)  <  1 . 5  VTir  uvwZ4, 
a  =  (D  -  R)/(V2M)  ,  D  =  V/H2+K2+L2 


Test  #3 


P  set  to  0  if  h  -  R  >  V5 TF2,  where  F2  =  4.262,  h  =  H,  K,  L, 
and  T  equals  the  corresponding  u,  v,  or  w. 


Test  #4 


P  set  to  1  if  P  =  P(R,H,K,L,M,M,M)  >  1  -  Z4 ,  where 

P  is  computed  from  Equation  10  or  Equation  11 


For  efficiency  it  is  important  to  try  to  reduce  the  number  of  times  Pc  is 
called  by  ELLCOV.  We  note  that  r  in  Equation  7  has  the  same  value  for  two  values 
of  z,  i.e.,  z  =  L  ±  (R  -  y),  where  0  <  y  <  R.  This  observation  in  many  cases 
reduces  by  half  the  number  of  times  Pc  needs  to  be  computed.  Indeed,  we  express 
Equation  6  as  the  sum  of  two  integrals,  the  first  has  the  limits  of  integration 
L  -  R  to  L  and  the  second  L  to  L  +  R.  In  the  second  integral,  the  substitution 
Z  =  2L  -  z  is  made  with  the  following  result 

1  fL 

P  =  -  |  [ E(z/w)  +  E[(2L  -  z)/w]  P  (r , H , K, u , v)  dz,  (23 

VSw  ->L-R 

or,  in  terms  of  Equation  8, 

i  rB5 

P  =  -  I  jexp(-t2)  +  exp[-(\/2L/w  -  t ) 2 ] j  P  (r,H,K,u,v)  dt,  (24 

J  A5  C 
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where 


R2  -  (L  -  \/2vt)2 


A5  =  (L-R)/(V2w) 


B5  =  L/(V2w). 


The  number  of  Gaussian  abscissae  used  to  integrate  Equation  24  is  fixed 
at  24,  and  in  some  cases  at  48.  Greater  accuracy  may  be  obtained  if  the  total 
integration  interval,  B5  -  A5  (=  R/(\/Y  w)),  can  be  reduced.  Indeed,  consider 
the  infinite  volume  Ve  bounded  by  the  planes  |z|  =  aw,  where  a  is  uniquely  defined 
by 


raw  f00  r°° 

'  (2TT\/2'rruvw)  l  E(z/w)  dz  l  E(y/v)  dy  I  E(x/u)  dx  =  1  -  e , 

J  —awj  J  —  oo  J  —oo 


£  >  0 , 


erfc  (a/VT)  =  e. 


The  entire  region  above  the  plane  z  =  aw  contains  e/2  of  the  distribution  F  given 
in  Equation  2.  Some  typical  values  for  a  as  a  function  of  e  are 


e  =  5 (-8 ) 


e  =  5/3(-9 ) 


a  s  3.8545\/2"  s  5.45109 
a  s  4.262VT  =  6.02738. 


In  ELLCOV  we  use  the  latter.  The  integration  over  the  sphere  S,  as  indicated  in 
Equation  6,  therefore  can  be  reduced  to  the  intersection  of  S  with  Ve  (SHVe).  In 
many  cases  this  permits  the  z  (or  t)  interval  of  integration  in  Equation  6  (or 
Equation  24)  to  be  shortened.  The  final  results  for  Equation  24  are  summarized  by 


L/V2  w 

a/\fl 


if  L  i  aw 


if  L  >  aw 


-a/Vz 

( L  -  R )  A/2  w 


if  L-R  <  -aw 


if  j  L  -  R  <  aw 


mi 
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0  if  L  <  aw 

.  .  (30 

1  otherwise. 

(Note:  P  set  to  0  if  L  -  R  >  aw. — Reduces  to  Test  #3,  a  = V?  F2) 

The  parameter  W9  is  used  in  ELLCOV  to  indicate  if  the  second  exponential 
in  Equation  24  can  be  ignored.  The  case  W9  =  1  implies  the  exponential  can  be 
dropped  and  W9  =  0  implies  it  must  be  retained.  The  case  for  W9  =  1  occurs  when 
the  integration  from  L  to  L  +  R  in  Equation  6  is  negligible  since  the  center  of 
S  lies  outside  S  D  Ve,  i.e.,  L  >  aw. 

A  further  effort  is  made  to  reduce  B5-A5  by  attempting  to  increase  A5. 

Let  the  integrand  in  Equation  24  be  denoted  by  G(t)  and  the  increasing  Gaussian 
abscissae  by  t  j ,  j  =  1,2,.. .,24.  If  G(A5)  <  1(-10)  then  G(tj)  is  evaluated  for 
increasing  j  until  G(tn)  >  1(-10).  Then  A5  is  set  to  tn-1  and  new  tj  are  gener¬ 
ated  based  on  the  new  A5. 

Since  G(t)  has  an  infinite  slope  at  t  =  (L  -  R)/(\/^2w),  it  may  happen  that 
accuracy  is  lost  when  G(A5)  and  G(tj)  differ  significantly.  Once  A5  is  deter¬ 
mined,  as  described  above,  if  Gft^)  >  1(~7)  a  24-point  Gaussian  integration  is 
carried  out  from  A5  to  t^ ,  followed  by  another  24-point  Gaussian  integration  frca 
tj  to  B5.  In  most  cases  GCt^)  <  l(-7),  in  which  case  only  the  24-point  Gaussian 
integration  from  A5  to  B5  is  used. 

Another  sensing  to  expedite  the  computation  determines  when  Pc  can  be  set 
to  one  from  some  t  to  B5.  Using  the  procedure  given  in  Reference  4  for  testing 
whether  Pc  >  1  -  e ,  we  have  Pc  >  1  -  1.3(-8)  if 

r  =  [R2  -  (L  -  vTwt)2]1/2  >  Vh2  +  K2  t  bs  e  G,  |  VTwt  -  L|  <  R,  (31) 

where 

s  =  max  (u,v),  b  =  6.02738. 

Therefore , 

|L  -  V2wt|  <  (R2  -  G2)1/?  =  U, 
or 

(L  -  U)/(V2w)  <  t  <  (L  +  U )  /( w ) ,  t  =  (L  -  U)/V2w.  (32) 


I  «,»  !,«  ».* 
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Again,  the  interval  of  integration  over  which  the  Gaussian  quadrature  is  needed 
can  be  reduced  when  Relation  31  holds  and  Equation  24  becomes 

P  =  l/VT  f  {  E(V2  t)  +  E[(2L  -  %/T  wt)/w]|  P  (r,H,K,u,v)  dt 
JAs 

r(L+U)/V2w 

+  i/a/7  (  E(V2t)  dt.  (; 

J (L-U) /y/l w 


We  also  note  that 


/•(L+U)  /\/2  w 

L/VT  l  exp(-t2)  dt  = 

-,(L-U)/V2w 


1/2  erfc [ (L  -  U)/\^Fw] 


if  (L  +  U )/y/l  w  >  F2 


|  l/2jerf[(L  +  U)/V^w]  -  erf[(L  -  U)/>/2w]j,  oth  erwise. 

Now  that  the  effective  range  of  integration  for  Equation  8  has  been  defined 
by  Equation  29  we  can  return  to  case  C.  If  A5  <  -  F2  then  it  is  possible  for  Pc 
to  remain  essentially  constant  over  the  effective  range  of  integration. 

Let  r(t)  =  yj  R2  -  (L  -  \/2wt)2.  Then  since  r(t)  is  an  increasing  function 
of  t,  we  take  Pc  as  a  constant  if 

|r(A5)  -  r(B5)|  <  l(-8)r(A5),  A5  <  -F2.  (-4.262) 

In  this  case  P  is  given  by 

P  =  erf (F2)  Pc|.5[r(A5)  +  r(B5) ] ,H,K,u,vj . 

After  much  experimentation  which  included  an  extensive  look  at  Simpson's 
rule,  the  24-point  Gaussian  integration  scheme  was  chosen  as  the  basic  quadrature 
procedure.  For  small  P,  a  lower  Gaussian  procedure  would  have  sufficed,  nc  er- 
theless  it  was  decided  for  simplicity  to  use  the  24-point  Gaussian  throughout. 

The  integrating  routine  is  called  GQUAD.  As  a  matter  of  interest,  GQUAD  allows 
for  the  interval  [A5,  B5]  to  be  subdivided  into  a  set  of  equal  subintervals  with 
the  same  Gaussian  order  formula  6,8,12,16,20,  or  24  applied  to  each  subinterval. 
This  feature  of  GQUAD  however  is  not  used  in  ELLCOV. 
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Extensive  testing  of  ELLCOV  showed  that  P  can  be  obtained  with  approximately 
8-decimal  digits,  even  though  we  conservatively  claim  6-decimal  digits.  There 
does  not  appear  to  be  any  significant  limitations,  for  realistic  values  of  the 
normalized  off-set  distance  V (H/u) 2  +  (K/v) 2  +  (L/w) 2,  and  the  quantities  u/v, 
v/w,  R/w  can  have  values  as  large  as  105. 

Below,  Table  1  contains  some  numerical  results  from  ELLCOV.  It  gives  values 
of  P  to  compare  with  approximations  given  previously  be  F.E.  Grubbs. a  listing 
of  ELLCOV  in  HP-BASIC  is  given  in  Appendix  C. 


R  =  1 


K  =  L  =  0 


a2  =  1/2 

a2  =  1 

a2  =  2 

a2  =  3 

u  =  v  =  w 

WH 

.890 

.607 

.313 

.204 

H  =  0 

P 

.888390 

.608375 

.317730 

.198748 

MC 

.885 

.615 

.320 

.210 

wm 

WH 

.807 

.530 

.307 

.205 

P 

.736060 

.501231 

.277180 

.179569 

■I 

MC 

.755 

.505 

.315 

.195 

WH 

.337 

.261 

.174 

.127 

H  =  1 

P 

.337133 

.269975 

.183329 

.  132299 

MC 

.353 

.320 

.205 

.135 

WH 

.000 

.001 

.004 

.006 

P 

3.0986(-3) 

1.59356(-2) 

3 . 3474 (-3) 

3.85359C-2) 

MC 

.000 

.020 

.035 

.040 

WH 

.000 

.000 

.000 

.000 

H  =  3 

P 

1.48C-7) 

7. 569(-5) 

1 . 7473 (-3) 

4 . 7661 (-3) 

MC 

.000 

.000 

.000 

.010 

u  =  2v  =  2w 

WH 

.868 

.630 

.388 

.279 

o 

II 

33 

P 

.876323 

.647542 

.375392 

.246952 

MC 

.860 

.670 

.400 

.250 

WH 

.739 

.554 

.364 

.269 

H  =  .5 

P 

.747518 

.569818 

.346921 

.233510 

MC 

.770 

.555 

.400 

.285 

WH 

.432 

.372 

.287 

.233 

H  =  1 

P 

.437088 

.385964 

.273711 

.  197402 

MC 

.435 

.405 

.290 

.235 

mm 

WH 

.019 

.055 

.092 

.104 

P 

3.08781 (-2) 

7 .57294 (-2) 

.  105474 

. 100705 

mi 

MC 

.040 

.090 

.  125 

.  105 

wm 

WH 

.000 

.001 

.011 

.023 

1 

P 

1 . 679 (-4) 

4.2614 (-3) 

2.11601 (-2) 

3 . 26805 (-2) 

■H 

MC 

.000 

.000 

.035 

.050 

a2  =  u2  +  v2  +  w2 


WH — Grubbs'  approximation,  MC — Monte  Carlo  (200  shots),  \ 

WH  and  MC  results  taken  from  [8],  [9].  P  computed  by  ELLCOV .  \ 
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TABLE  1  (Cont) 


R  =  1 


a2  =  1/2 


.878 

.878302 

.890 


.738 

.735871 

.735 


.378 

.378261 

.380 


.005 

9 . 5532  (-3) 
.015 


.000 

4.86(~6) 

.000 


.740 

.869274 

.774 


.739 

.747340 

.750 


.447 

.452305 

.452 


WH  .000 

P  3.428(-4) 

MC  .001 


a2  =  1 


.618 

.625938 

.635 


.524 

.531421 

.520 


.308 

.319040 

.355 


.021 

3 . 40282 (-2) 
.040 


.000 

5 . 545 (-4) 
.000 


u  =  2v  =  3w  | 

.640 

.659413 

.662 


.568 

.586033 

.584 


.392 

.409374 

.408 


.025  .068 

3.96655(-2)  9.20037(-2) 

.044  .087 


.002 

6 . 6584  (-3) 
.009 


a2  =  u2  +  v2  +  w2 


L  =  0 

o2  =  2 

o2  =  3 

.350 

.351727 

.365 

.236 

.229160 

.275 

.314 

.315323 

.395 

.219 

.211615 

.270 

.221 

.226834 

.265 

.  171 
.  166578 
.210 

.044 

5.94781 (-2) 
.050 

.055 

6 . 36591 (-2) 
.070 

.002 

6 . 0278 (-3) 
.000 

.006 

1 . 2624 1 (-2) 
.010 

.415 

.407256 

.412 

.318 

.279152 

.275 

.390 

.378893 

.366 

.299 

.265242 

.266 

.312 

.305022 

.306 

.259 

.227522 

.223 

.  109 
. 127553 
.122 

.123 
.  123078 
.  120 

.015 

2 . 94238 (-2) 
.029 

.030 

4 . 40764 (-2) 
.046 

WH — Grubbs'  approximation,  MC — Monte  Carlo  (200  shots), 

Last  MC  Group  (u  =  2v  =  3w)  1000  shots, 

WH  and  MC  results  taken  from  [8],  [9].  P  computed  by  ELLCOV. 
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III.  COMPUTATION  OF  R  (The  inverse  problem) 

In  this  section  we  describe  a  procedure  for  finding  R  given  P,  H,  K,  L,  u, 
v,  w.  For  fixed  H,  K,  L,  u,  v,  w  the  radius  R  is  a  monotone  increasing  function 
of  P  and  thus  unique  for  a  given  P.  The  procedure  for  finding  R  requires  the 
computation  of  P,  as  described  in  the  previous  section.  Since  that  basically 
involves  a  time-consuming  numerical  triple  integration  a  large  effort  is  made 
to  obtain  good  early  approximations  for  R.  Once  such  approximations  are  found 
uniform  stepping,  halving,  Regula-f alsi ,  and  King's  root  finding  procedure**  are 
used  to  refine  the  estimates  for  R.  The  objective  is  to  obtain  R  correctly  to 
approximately  6-decimal  digits. 

In  the  discourse  some  details  are  omitted.  They  may  be  obtained  from  a 
detailed  study  of  INVELLCOV,  the  BASIC  program  for  computing  R.  A  listing  is 
given  in  Appendix  C  and  a  flowchart  is  shown  in  Figure  2  on  page  18. 

The  unknown  value  of  R  will  always  be  contained  in  a  known  interval. 

Initially  crude  lower  and  upper  bounds,  Rmin  and  Rmax,  are  found  for  R. 

Let  I  H  max[3  W/2  Puvw,  ( V  tt/2  PM)3],  M  =  max(u,  v,  w).  Then 

I1/3  if  P  i  1/2 

R  £  R  .  =  I1/3  if  P  >  1/2  and  I  >  D3  =  (H2  +  K2  +  L2)3/2  ( 34 ) 

min 

D  if  P  >  1/2  and  I  <  D3. 

The  derivation  of  Equation  34  is  similar  to  the  methods  used  for  deriving  Tests  #1 
and  t 2  of  Section  II.  See  Appendix  B. 
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For  Rmflx,  two  arrays  A6(j),  B6(j),  j  =  1,2, ...10  are  used. 


A6(l )  =  5 (-6) 
A6(2)  =  l(-4) 
A6(3)  =  .01 
A6(4 )  =  .1 
A6(5)  =  .3 
A6(6)  =  .6 
A6 (7 )  =  .9 
A6(8 )  =  .999 
A6(9)  =  .999999 
A6(10)  =  1 


B6(l )  =  .0266 
B6(2)  =  .0723 
B6(3)  =  .339 
B6(4)  =  .765 
B6 ( 5 )  =  1.1933 
B6(6 )  =  1.717 
B6 (7 )  =  2.5005 
B6 ( 8 )  =  4.0335 
B6 (9 )  =  5.538 
B6(10)  =  10 


The  element  B6(j)M  gives  the  radius  of  a  sphere  centered  at  (0,0,0)  for  which 
A6( j)  =  P(B6( j)M, 0,0,0, M,M,M).  Then 


R  <  R  =  D  +  B6(J)M, 
max 


M  =  max(u,v,w) , 


where  J  is  the  minimum  integer  j  for  which  P  i  A6(j).  The  plausibility  of  Equa¬ 
tion  35  is  easily  seen  from  its  2-dimensional  analog.  In  Figure  1,  the  inner 
ellipse  contains  P  of  the  distribution.  The  outer  ellipse  with  semi-axes  qu,  qv 
contains  A6(J)  of  the  distribution.  Then,  from  the  figure,  one  easily  concludes 
R  <  V  H2  +  K2  +  qs,  where  s  =  max(u,v)  and  q  corresponds  to  B6(J). 


FIGURE  1:  USED  WITH  (35) 
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^min  and  Rmax  es  tablish  crude  bounds  for  R.  An  improvement  is  generally 
obtained  by  using  an  estimate  for  R,  call  it  Rg,  given  by  F.E.  Grubbs.8*^  His 
estimate  depends  on  a  percentage  point  of  the  chi-squared  distribution  which  is 
available  through  the  subprogram  GAMINV.  This  subprogram  is  described  in  Refer¬ 
ence  7.  The  quantity  actually  given  by  Grubbs  is  Rg2  which  estimates  R2  and  is 
given  by 


Rg2  =  W4{[x(4/V5,P)  -  4/V5]W2  +  Nj, 


where 


W4  =  u2  +  v2  +  w2 
N  =  1  +  D2/W4 


•■-■lit'-®]*  5H»1*i  ['*-©11 
15  ['  *’®1*  5  [■  *’®‘j  *ii‘  *>|»|’]l 


W2  =  T5/(2V4) 

V5  =  T52/V43, 

and  x(A,P)  satisfies 


(37.1) 

(37.2) 


(37.3) 


(37.4) 


(37.5) 

(37.6) 


u; 


'X(A.P) 


exp(-t)  t  dt,  A  =  4/V5 .  [1;  p.  260] 


Rmin2  <  RR2 


<  Rg  <  Rn 


then  Rg  is  accepted  as  an  improved  estimate  for  R.  If  Equation  39  holds  and 
P(Rg)  <  P,  then  Rg  is  stored  in  Q4  and  Rmax  is  stored  in  Q5;  if  P(Rg)  >  P,  then 
Rg  is  stored  in  Q5  and  Rm^n  is  stored  in  Q4.  Generally  Q4  and  Q5  bound  R  with 
Q4  <  R  <  Q5  (When  King's  procedure  is  used  Q4  and  Q5  still  bound  R,  but  their 
roles  as  lower  and  upper  bounds  may  be  reversed.).  Refinements  to  Q4  and  Q5  are 
obtained  by  using  the  procedures  mentioned  above.  Some  of  the  details  are  shown 
in  the  flow  chart  of  Figure  2.  A  flow  chart  for  King's  root  finding  method  is 
given  in  Figure  3  on  page  19. 
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NOTATION  FOR  FLOWCHART  OF  FIGURE  2  AND  BASIC  LISTING 


Rt  denotes  true  value  of  R  for  a  given  value  P  =  P3.  P3  =  P(Rt) 

C2  =  10  X8  (1  -  P3) ,  P3  <  .5;  C2  =  10  X8  P3,  P3  >  .5 

E4  =  0,  1 (-6 )  <  P3  <  .999999;  E4  =  1,  l<-6)  >  P3  >  .999999 

E5  =  0,  Less  than  30  iterations  to  find  R3;  E5  =  1 ,  30  iterations  are  used  to 

find  best  estimate  for  Rt.  Final  estimate  R3  may  not  be  accurate  to 
6-decimal  digits. 

Hx  -  H  of  text;  Hy  -  K  of  text;  Hz  -  L  of  text. 

16  -  R3  with  PA  <  0 

17  -  R3  with  PA  >  0 

18  -  P(R3),  latest  value  of  P 

19  -  Increment  used  in  stepping  procedure.  19  =  ±.2  (Q5  -  Qi+). 

K0  =  H2  +  K2  +  L2;  K2  =  VK0 

P  -  Previous  value  of  18  -  P3 
P3  -  Input  value  of  P 

PA  -  P(R3)  or  P(R3)  -  P3 

P5  -  P(Q4)  -  P3 


-  P(Q4)  -  P3 

-  P(Q5 )  -  P3 

-  P5  or  P6 

-  Latest  value  of  P  >  P3 

-  Latest  value  of  P  <  P3 

-  See  Section  III.  Latest  estimate  for  Rt  such  that  P(Qi+)  -  P3  <  0. 

-  See  Section  III.  Latest  estimate  for  Rt  such  that  P(Qs>  -  P3  >  0, 

-  Contains  number  of  iterations  used  to  estimate  Rt. 

-  Latest  estimate  for  Rt 

-  See  Section  III.  Estimate  from  Grubb's  approximation  for  Rt . 

-  See  Section  III.  Lower  estimate  for  Rt  from  ELLRC. 

-  The  latest  estimate  for  Rt. 

-  The  next  to  latest  estimate  for  Rt . 

-  Starting  lower  bound  for  Rt.  Rmax  -  Starting  upper  bound  for  Rt  . 

-  u  of  text 

-  v  of  text 

-  w  of  text 

-  max(u,  v,  w) 

-  Number  of  iterations  used  to  obtain  Rc  using  ELLRC. 

-  See  Section  III.  Used  to  exit  INVELLCOV  if  | PA (R3 )  -  P3  |  <  X8. 
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Input:  xQ,  fQt  fj,  (fQf1  <  0) 


Replace 

(x0,f0)by(x1,fi), 

(Xpf  j)by(x2,f2)  . 


START 


Use  (*)  for 
new  X2 
Compute  f2 


I  ntercliange 
(x0,f  JantKxj.fj) 


flf2<0? 


Use  ( * )  to 
get  new  x2 
Compute  f? 


x_  =  latest  estimate  for  x,  f  =  f (x  ) ,  f(x) 
l  n  n 


FIGURE  3.  KING'S  ROOT  FINDING  METHOD  FOR  f(x) 
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In  most  cases  Rg  gives  a  very  good  estimate  for  R.  In  fact,  Rg  in  many  cases 
can  be  accepted  as  the  final  estimate  for  R,  with  no  iterations  required.  If 
however  Inequality  39  does  not  hold,  (the  right  hand  side  of  Equation  36  may  even 
be  negative)  then  we  proceed  as  follows. 

The  constant  Rp  is  determined  so  that 

fL+R  fK+R  rH+R 

P  -  l  I  I  F(x,v,z,u,v,w)  dx  dy  dz.  (40) 

J L-Rp  K-Rp J H-Rp 

The  subroutine  ELLRC,  using  the  Newton-Raphson  method,  is  used  to  compute  Rp. 

The  right  hand  side  of  Equation  40  yields  the  probability  P  over  the  smallest 

cube,  with  center  at  (H,K,L)  and  sides  parallel  to  the  x,  y,  and  z  axes,  which 
contains  the  sphere  of  radius  Rp.  Hence  Rp  <  R  <  \/3  Rp;  if  Rp  >  Rmjn,  then 

Q4  =  Rp,  and  if  y/~3  Rp  Rmax  then  Q5  =  \/3  Rp. 

In  ELLRC  the  maximum  number  of  iterations  for  Rp  is  set  at  40.  A  maximum  of 
37  iterations  has  been  observed. 

Again,  once  Rp  is  computed,  refinements  are  made  by  using  halving,  stepping, 
Regula-falsi  and  King's  method  to  obtain  a  final  estimate  for  R. 

The  iterations  for  R  are  generally  terminated  when  Q4  and  Q5  agree  to  at 

least  6  significant  digits  or  if  the  latest  estimate  for  R,  say  R3 ,  yields  P(Rq) 

such  that  |  P  -  P(R()[  £  Xg,  where  Xg  is  a  small  number  which  depends  on  P.  For 
example,  if  P  <  l(-5)  then  Xq  =  5(-10),  whereas  if  .001  <  P  "  .999995  then 
X p  =  5(-7).  When  exit  conditions,  some  of  which  have  just  been  described,  are 
satisfied,  then  one  more  application  of  Regula-falsi  is  carried  out  before  exiting 
INVELLC0V.  See  Figure  2  for  more  details,  page  18. 

The  maximum  number  of  iterations  for  R  in  INVELLC0V  is  set  at  30.  Extensive 
checking  has  never  shown  more  than  20  iterations  for  any  case.  The  average  number 
of  iterations  is  between  4  and  5.  When  P  is  very  small  (<.001)  more  than  the 
average  number  of  iterations  is  generally  needed. 

INVEI.LCOV  can  be  used  with  any  reasonable  input  values  of  H,  K,  L,  u,  v,  w. 
The  input  variable  P  should  satisfy  jp  -  1/2 j  <  .499999.  If  P  is  in  (0,1)  and 
does  not  satisfy  the  above  inequality,  INVELLC0V  may  still  give  a  result  for  R, 
but  there  is  no  assurance  near  6  digit  accuracy  will  be  obtained. 

Table  2  below  contains  some  numerical  results  from  INVELLCOV.  It  indicates 
the  wide  range  of  applicability  of  INVELLCOV. 
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TABLE  2 


Item 

u/v 

v/w 

R/w 

H/u 

K/v 

L/w 

p 

1 

1/2 

2/3 

3.183274/6 

5/2 

5/2 

10/3 

5  (-6) 

2 

1/2 

2/3 

10.28906/6 

5/2 

5/2 

10/3 

5  (-3) 

3 

1/2 

2/3 

16.60908/6 

5/2 

5/2 

10/3 

.  1 

4 

1/2 

2/3 

23.39931/6 

5/2 

5/2 

10/3 

.5 

5 

1/2 

2/3 

30.48622/6 

5/2 

5/2 

10/3 

.9 

6 

1/2 

2/3 

40.77124/6 

5/2 

5/2 

10/3 

.999 

7 

1/2 

2/3 

48.44240/6 

5/2 

5/2 

10/3 

.999995 

8 

l(-3) 

.05 

l.lKlll 

4000 

20 

5 

5  (-6) 

9 

l(-3) 

.05 

2.629994 

4000 

20 

5 

5(-3) 

10 

l(-3) 

.05 

3.855994 

4000 

20 

5 

.  1 

11 

l(-3) 

.05 

5. 103200 

4000 

20 

5 

.5 

12 

l(-3) 

.05 

6.364060 

4000 

20 

5 

.9 

13 

l(-3) 

.05 

8.155921 

4000 

20 

5 

.999 

14 

l(-3) 

.05 

9.472573 

4000 

20 

5 

.999995 

15 

2/5 

5 

6.184626 

0 

1 

10 

5  (-6) 

16 

2/5 

5 

8.185500 

0 

1 

10 

5  (-3) 

17 

2/5 

5 

9.733491 

0 

1 

10 

.  1 

18 

2/5 

5 

11.72659 

0 

1 

10 

.5 

19 

2/5 

5 

15.45659 

0 

1 

10 

.9 

20 

2/5 

5 

22.95502 

0 

1 

10 

.999 

21 

2/5 

5 

29.02551 

0 

1 

10 

.999995 

22 

K-5) 

2(9) 

2. 1376073(4) 

0 

l(-4) 

0 

5  (-6) 

23 

l(-5) 

2(9) 

1.2533239(7) 

0 

l(-4) 

0 

5  (-3) 

24 

l(-5) 

2(9) 

2.5132270(8) 

0 

K-4) 

0 

.  1 

25 

l(-5) 

2(9) 

1.3489795(9) 

0 

K-4) 

0 

.5 

26 

l(-5) 

2(9) 

3.2897073(9) 

0 

K-4) 

0 

.9 

27 

K-5) 

2(9) 

6.5810535(9) 

0 

K-4) 

0 

.999 

28 

K-5) 

2(9) 

9.1295755(9) 

0 

K-4) 

0 

.999995 

Values  of  R  as  computed  by  INVELLCOV. 
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Derivation  of  Equations  10,  11  and  13-17 

Equations  10  and  11  refer  to  Case  A  in  Section  II.  We  have 
u  =  v  =  w|  D  =  VH2  +  K2  +  L2, 


P  =  J*F(x,y,z,u,v,w)  dx  dy  dz, 


where  the  integration  is  carried  out  over  the  interior  of  the  sphere 
S  :  (x  -  H)2  +  (y  -  K)2  +  (z  -  L)2  =  R2 . 

Let 

x  =  H  +  t  cos  9  sin  <f> 

y  =  K  +  x  sin  0  sin  4> 

z  =  L  +  x  cos  <f> , 

O<t£R,O£.0£  2tt,  0  l  <f>  £  7i. 

Then  Equation  A-2  becomes 

i  fR  cv  c2v 

p  -  - 1  I  1  E(D/u)  exp/-  —  [H  cos  0  cos  <J>  +  K  sin  0  sin 

2tt  V'27uvwJ0Jc!>=0  •'0=0  <  u2 

+  L  cos  $  ]  |  E(t/u)  t2  sin  <f)  d0  d<j>  dx . 

Taking  advantage  of  spherical  symmetry,  we  take  H  =  K  =  0,  L  =  D.  Hence 

i  fR  c 

P  =  — - I  1  E(D/u)  exp(-Dx  cos  4>/u2)  E(x/u)  sin  4>  d<}>  x2  dx  . 

Vl tT  u3  '  0  ^0 

With  cos  <J>  =  s , 


(A—  1 ) 


(A-2) 


(A-3) 


E(D/u)  exp(-Dx  cos  4>/u2)  E(x/u)  sin  4>  d<}>  x2  dx . 


(A-4) 


fR 

r1 

l  E (D/u)  E(x/u) 

l  exp(-Dxs/u2)  ds 

■'O 

-J-l  J 

where 


n/2tT  u3  *  0 


exp(-Dxs/u2)  ds  =  [exp(Dx/u2)  -  exp(-Dx/u2) 


m 

I 

i 


! 


I 


i 

Emu 
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Thus 


P  = 


i  rR 

-  I  (e[(D  -  t)/u]  -  E [ (D  +  t ) /u ] |  t  dx, 

V2k  Du  ^0 


and 


■  i  -  -(s')  •  *  5  ft**)  -  -f-^i  ^ 


(A-5) 


If  D  =  0,  then 


P  =  erf 


/-M- V1  s 

\V2u)  y  77  u 


-  E(R/u) ,  D  =  0. 


(A-6) 


Now  let 


Y  = 


2DR/u2 


2DR/u" 


Using  the  result  in  Equation  A-5  gives  Equation  10;  then  Equation  A-6  or  Equa¬ 
tion  11  follows  since  Y  =  E(R/u)  when  D  =  0. 


Derivation  of  Equations  13-17 


We  have 


u  =  v  and  H  =  K  =  0 


u  ^  w. 


(A-7) 
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Therefore 


P  2  [erf  L4  erf  L5]  T>  L4  ^  w  ’  L5  "  ^  w 


e(Vr2  -  L2/u)  fL+R 


E(z/w )  exp£-  -j(2Lz  -  z2)/u2Jdz 


(A-9) 


(A- 10) 


Two  cases  are  now  treated  separately,  namely  u  >  w  and  u  <  w. 

Let  u  >  w.  Then  Equation  A-10  becomes,  after  straight-forward  algebra, 


F1=[a(L+R)+B]/\/2 


T  =  -  exp 

y/Ti 7  w 


• -n© 


=  [a(L-r)+B  ]/y/2 


exp(-t2)  dt,  (A— 11) 


where 


a  =  L2/w, 


We  also  have 


l2  =  V  I  1  “  (w/u)2 


(L  +  R)L  , 

F.  2  [a(L  +  R)  +  8]  /  \fl  =  - -  +  -  =  -  (L/L.  +  RL  )  (A-12.1) 

yfl  w  VT  L2u2  V2  w 


Sn  =  [  a(L  -  R)  +  B  ]  /  y/l  =  — - —  (L/L.  -  RL.),  (A-12.2) 

0  V2w  2  2 

with  L?  given  by  Equation  19.1.  Then  from  Equation  A-9,  using  Equation  A-ll 
and  Equation  A- 12,  we  get 

P  =  -j  jerf  L^  -  erf  L^  -  E(R/u)  E[L/(uL2)]  (erf  F1  -  erf  SQ)|  (A-13) 

which  is  Equation  13. 

Another  form  for  Equation  A-13,  namely  Equation  14,  is  obtained  by  replacing 
erf  Fj  -  erf  by  erfc  Sq  -  erfc  F^,  and  then  multiplying  erfc  Sq  by 
expCS^2)  expC-S^2),  and  erfc  F^  by  exp(F^2)  exp(-F^2).  Combining  the  results 
with  the  remaining  exponentials  in  Equation  A-13  gives  Equation  14,  where  use 
is  made  of  the  relation 

y  | (R/u)2  -  L2 /  (u2  -  w2)  +  [a (L  ±  R)  +  B]2j  =  y[  (L  ±  R)/w]2. 


P 

s$ 

m 

»»>« 


i 


W 

JlirJi 


NSWC  TR  87-27 


When  L  =  0,  Equation  A-13  reduces  to 


P  =  erf(R/VT  w)  -  (1/ L^)  E(R/u)  erf  (RI^/V^w) . 

Now  consider  w  >  u.  Starting  with  Equation  A-IO,  we  have 

_  1  R  |- .  -i 

T  =  E(R/u)  E(L/Vw2  -  u2)  -  1  exp  I  y  (az  -  8)2  dz. 

■v/2tT  w  •'L-R  ^  ^ 


(A-14) 


(A-15) 


Note  here  that  the  argument  of  the  exponential  in  the  integrand  is  positive. 
Proceeding  as  above,  T  is  obtained  in  terms  of  Dawson's  integral  denoted  here 
by  daws (x) . 


Vir  L5 


E(R/u)  e|l/  \/(w2  -  u2  )  J  Jexp  (F^  2  )daws  (F^  )  -  exp  (Sq 2  )daws  (Sq  , 


daws(x)  =  exp(-x2)  [  exp(t2)  dt,  daws(-x)  =  -daws(x),  [1;  p.  298],  [3] 


We  also  have  using  Equation  A- 12 


[R(w2  -  u2)  -  Lu2] 


y/2  uw  V  I  u2  -  w2 


y/2  uw  v  |  u2  -  w2  I 


[R(w2  -  u2)  +  Lu2 ] 


Thus  using  the  relation 


1  |(-)2  +  - — - - -  [ R2  (w2  -  u2)2  ±  2RLu2  (w2  -  u2)  +  u4L2]l 

2  (w2  -  u2)  u2w2(w2-u2)  } 


- 1|(!) 


2(l  _  (*2  -  u2)\  ,  2RL  +  L2 
\  w2  /  w2  (w2  -  u2 ) 


1  -  Hi 


■[— T 

L  vTwJ 


Equation  A-15  reduces  to 


T  =  - —  e(— - — ) [daws (F,  )  -  exp(-  2RL/w2)  daws(Sn)], 

2L2  ^  \  w  1  0 

which  agrees  with  the  corresponding  term  in  Equation  16. 


(A— 16) 


A-6 
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For  L  =  0,  with  w  >  u,  it  is  easy  to  conclude  from  Equation  A-9  and 
Equation  A- 16  that 

0  1 

P  =  erf  L. - - —  E(R/w)  daws (RL./V2 w)  , 

4  /—  l 

Vir  2 

which  is  Equation  17. 


jffl 

teJ:: 
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Discussion  of  Tests  #1,  // 2,  #3,  #4 


Test  // 1 :  From  Equation  4  with  E(x/u)  E(y/v)  E(z/w)  <  1 ,  we  have 


P  <  R3/(2ttV2tT  u  v  w)  <  Z4 ,  Z4  =  2.5(-8), 


R3  <  1.5  uvw  Z4  (=>  P  <  Z4). 


(  B  —  1 ) 


The  relations  in  Inequality  B-l  lead  to  Test  #1. 


Also,  from  Equation  6,  setting  E(z/w)  P^  =  1  with  w  =  max(u,v,w) 


P  <  —— —  [(L  +  R)  -  (L  -  R)] 
V2tt  M 


<  Z4,  M  =  max(u,v,w) , 


R  <  y-  Z4M. 


(B-2) 


The  last  inequality  and  Inequality  B-l  are  used  in  Equation  34. 


Test  #2:  Note  if  D  >  R  then  P  <  1/2  and  the  sphere  S  does  not  contain  the 
origin.  Therefore,  we  have  for  D  >  R 


E (x/u)  E(y/v)  E(z/w)  <  E(x/M)  E(y/M)  E(z/M)  <  E[ (D  -  R)/M] 


and  it  follows 


P  <  - 1 -  f  E 

2tt  V2tt  uvw  ’'S 


D  -  R 
M 


dx  dy  dz  =  - -  ■ —  e(~— i  Z4 .  (B-3) 

6it  \/ 2tt  uvw  ' 


Then  Test  #2  follows  directly  from  Inequality  B-3. 


Test  #3:  This  follows  by  recognizing  that  the  sphere  S  is  located  outside 
the  region  which  contains  1  -  Z4  of  the  distribution.  It  is  shown  by 
Equation  26  and  Equation  27  that  if  a  is  determined  by  erfc(a/-\/2)  =  e  then 
the  entire  region  0  above  the  plane  z  =  a  w  contains  c/2  of  the  distribution. 
Since  z  =  L  -  R  >  0  locates  the  point  of  the  sphere  closest  to  the  xy-plane, 
and  if  L  -  R  >  y/l  w  F2(Vr2  F2  =  a),  then  the  sphere  S  must  be  entirely  contained 


in  the  region  U.  Here  L  has  been  used  for  h,  and  w  has  been  used  for  T. 


W 

& 

i 


m 

M 

m 


1 


05 


IS 


■  **>  *>  *«*  v  v  v*  "v  ^  "v  -V<  .v*  h 


kl 


NSWC  TR  87-27 


Test  #4 :  This  test  is  based  on  the  fact  that  if  P  ~  1  and  u,  v  or  w  increases, 
then  for  fixed  R,  P  must  decrease.  Thus  if  we  replace  u,  v,  and  w  by  M  to  obtain 


P,  then  P  <  P.  If  P  >  1  -  Z4,  then  P  >  1  -  Z4. 


A  heuristic  argument  follows  that  indicates  if 


R2  >  H2  +  K2  +  L2, 


P  =  P(R,H,K,L,u,v,w)  >  P(R,H,K,L,M,M,M,),  M  =  max(u,v,w) 


Without  loss  of  generality,  assume  w  <  M  and  let 


P  =  P0  =  P(Rq ,H,K,0,u,v,w) . 


Then  it  is  easy  to  show  from  Equation  8  (with  L  =  0)  that 


3P0  ,8/  V2  w  3PC  8f 

T“  '  J-R/V2wlr3" 


dt  <  0, 


where 


3P  (C,H,K,u,v) 


'R2  -  2w‘  V  and 


Thus  from  Equation  B-7 , 


P  =  P0(R0,H,K,0,u,v,w)  >  P(R0,H,K,0,u,v,M) 


Then  by  continuity,  we  have  for  sufficiently  small  6 


P  =  P(R,H,K,S ,u,v,w)  >  P(R,H ,K,6 ,u ,v,M) . 


Now  consider  a  sequence  of  spheres  js^j  with  centers  (H,K,5j)  and 
radii  R^,  such  that  6i  is  an  arbitrary  increasing  set  of  positive 
numbers  with  0  <  5^  <  L,  R^  i  R^i  R  and 


P  =  P(Ri,H,K,6i,u,v,w). 


Note  that  by  Inequality  B-4  all  contain  the  origin. 


(B— 4 ) 


(B— 5  ) 


(B-6) 


(B-7) 


(B-8) 


(B-9) 


(B- 10) 


(B- II) 
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If  P  =  P(R^,H,K,6^,u,v,w)  >  P(R^,H,K,6^,u,v,M)  holds  for  all  elements  S^, 

then  a  simple  argument  establishes  Inequality  B-5.  Suppose  however  for  some 

R.  =  R,  6.  =6  the  last  inequality  does  not  hold.  Since  Inequality  B-9 
^  ^  A  & 

always  holds,  there  must  exist  by  continuity  a  sphere  S  with  R,  <5, 

_  *  a  -  *  * 

(R  <  R  <  R  <  R,  6  <  6  <  6  <  L)  such  that 

P  =  P(R,H,K,*,u,v,w)  =  P(R,H,K,*,u,v,M) .  (B-12) 

But  Equation  B-12  is  impossible,  since  it  implies  that  the  two  different 

cumulative  normal  distributions  over  the  same  sphere  can  give  the  same  value 

of  P.  Indeed,  denote  the  two  different  density  functions  of  Equation  (B-12) 

by  I  and  T  .  These  functions,  with  3  independent  variables  and  hence  not 
w  M 

easily  visualized,  intersect  along  a  surface  T  bounded  by  a  simple  closed 

curve.  From  Equation  B-12  and  Inequaltiy  B-4  with  the  integrations  of  I 

a  ^ 

and  I..  carried  out  over  the  interior  of  S,  we  have 
M 

P  ’  +  Av  *  JM  +  V 

where  J  (Jw)  denotes  the  integral  of  I  (I  )  over  that  part  of  S  which  contains 
w  M  w  M 

the  origin  and  is  bounded  by  T;  A^CA^)  denotes  the  remaining  contribution  from 
the  integration  of  I  (1^)  over  that  part  of  S  bounded  by  T  and  not  containing 
the  origin.  For  w  <  M,  then  I  >  I  at  the  origin  and 

W  M 

J  >  J.„,  or  J  =  *1  +  a  a  >  0;  (B— 14) 

w  M  w  M 

Aw  "  V  0r  AW  =  ^  -  b  b  -  °-  (B_15) 

But  a  >  b,  since  the  integration  of  I  and  1^  over  all  xyz-space  must  equal 

one.  This  leads  to  a  contradiction  of  Equation  B-13. 

This  argument  is  easily  seen  in  the  one  dimensional  case.  We  use  the 
same  notation  as  above  for  the  corresponding  quantities  in  one  dimension. 

See  the  figure  below. 
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FIGURE  4:  NORMAL  CURVES  WITH  UNEQUAL  STANDARD 
DEVIATIONS 


D(t/y)  -  -  E(t/y) 

•n/2tT  v 


r+)  =  f  +  D(t/w)  dt,  jm(c,t+)  =  f  + 
D(t/w)  dt,  =  J  D(t/M)  dt. 


Dft/M'*  dt 


D(t/M)  dt. 


Clearly  with  c  <  0,  d  >  0  (corresponding  to  +  L^)  we  have 

J  >  J  ,  A  <  A  ,  and  J  (c,d)  >  J  (c,d). 
w  M  w  M  w  M 

Thus  by  the  contradiction  above  we  conclude  that 
P(R,H,K,L,u,v,w)  >  P(R,H,K,L,u,v,M) 

Then  Inequality  B-5  easily  follows  by  carrying  out  the  same  argument  using  the 
min(u,v),  say  v,  in  place  of  w  to  show  the  last  inequality  below, 

P(R,H,K,L,u,v,w)  >  P(R,H,K,L,M,v,M)  >  P(R,H,K,L,M,M,M) . 
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APPENDIX  C 

H-P  BASIC  LISTINGS  OF  ELLCOV  AND  INVELLCOV 
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5  !  ARRAYS  X<*),  Y  <  *  )  ,  B  <  *  )  MUST  BE  CALLED  WITH  BASIC  SUBROUTINE 

6  !  IN  MAIN  PROGRAM  IN  ORDER  TO  USE  ELLCOV  OR  INVELLCOV. 

0  DIM  S2  <  20  >  ,  Hx<  20  >  ,  H  <  3 ) ,  S  <  3  ) ,  A5  <  3  )  ,  B5  <  3  ) ,  T5  <  3 )  ,  P  <  7  ) ,  W8  <  3  ) 

5  COM  X<43>, Y<43>, A6< 10) ,B6< 10) ,Q9, V,A7,P4, Ve,T6,Pl,T9 
0  !  GAUSSIAN  CONSTANTS  STORED  IN  "GAUSS".  USED  IN  "SUBGquad" . 

5  X< 1 )  =  . 238619186083 
0  X<2)=. 661209386466 

5  X<3)=. 932469514203 

0  X<4>=. 183434642496 

5  X<5)  =  . 525532409916 

0  X(6)=. 796666477414 

5  X<7)=. 960289856498 

0  X(8)=. 12523340851 1 

5  X<9)=. 367831498998 

0  X< 10)=. 587317954287 
5  X< 1 1 )=. 769902674194 
0  X< 12)=. 9041 1725637 
5  X< 13)=. 981560634247 
0  X<14)=9. 50125098376E-2 

5  X(15)v 281603550779 

0  X< 16)=. 458016777657 
5  X< 17)=. 617876244403 
0  X< 18)=. 755404408355 
5  X( 19)=. 865631202388 
0  X<20)=. 944575023073 

5  X(21 )=. 989400934992 

0  X<22)=7. 6526521 1335E-2 

5  X(23)=. 227785851 142 

0  X(24)=. 373706088715 

5  X<25)=. 510867001951 

0  X<26)  =  . 636053680727 

5  X(27)=. 74633190646 

0  X<28)=. 8391 16971822 

5  X<29)=. 912234428251 

0  X(30)=. 963971927278 

5  XC31 )=. 993128599185 
0  X  (  3  2 )  =  6 . 405689286 2 6E-2 
5  X<33)=. 191 1 18867474 

0  X<34)  =  .  315042679696 

5  X < 35 )=. 433793507626 
0  X<36)=. 545421471389 

5  X<37)=. 648093651937 

0  X(38)=. 740124191579 

5  X(39)=. 820001985974 

0  X(40)=. 886415527004 

5  X(41 )=. 938274552003 

0  X<42)=. 974728555971 

5  X<43)=. 995187219997 

0  Y(1 )=. 467913934573 

5  Y<2)=. 360761573048 

0  Y<3)=. 171324492379 

5  Y<4 )=. 362683783373 

0  Y<3>=. 313706645878 

5  Y(6)=. 222381034453 

0  Y(7)=. 10122853629 

5  Y<8)=. 249147045813 

0  Y(9)=. 233492536538 
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2205  Y< 10)=. 203167426723 
2210  Y< 1 1 )=. 160078328543 
2215  Y<12)=. 106939325995 

2220  Y< 13)=4. 71753363865E-2 
2225  Y< 14)=. 189450610455 
2230  Y(15)=. 182603415045 

2235  Y<16)=. 169156519395 

2240  Y< 17)=. 149595988817 
2245  Y< 18)=. 124628971256 
2250  Y< 19)=9. 515851 16825E-2 
2255  Y ( 20 ) =6 . 22535239386E-2 
2260  Y(21 )=2. 715245941 18E-2 

2265  Y<22)=. 152753387131 

2270  Y(23)=. 149172986473 

2275  Y<24)=. 142096109318 

2280  Y<25)=. 131688638449 

2285  Y<26)=. 1 18194531962 

2290  Y(27)=. 1019301 19817 

2295  Y ( 28 ) =8 . 32767415767E-2 
2300  Y<29)=6. 26720483341E-2 

2305  Y<30)=4. 06014298004E-2 

2310  Y<31 )=1 . 76140071392E- 2 

2315  Y(32)=. 127938195347 

2320  Y<33)=. 125837456347 

2325  Y(34)=. 121670472928 

2330  Y<35)=. 1 15505668054 

2335  Y<36)=. 1074442701 16 

2340  Y( 37 ) =9 . 76186521041E-2 
2345  YC38J=. 086190161532 

2350  Y<39)=7. 334648141  1  IE-2 

2355  Y<40)=5. 92985349154E-2 

2360  Y(41 )=4. 42774388174E-2 

2365  Y(42)=2. 853 1 3886289E-2 

2370  YC43)=. 0123412298 


23"’5 

86(1 )=5E-6 

2380 

86  C  2 )  =  . 000 1 

2  3  85 

86  <  3 )  =  . 0 1 

2390 

86(4)=. 1 

2395 

86  <  5 )  =  . 3 

2400 

86(6)=. 6 

2405 

86(7)  =  .  9 

2410 

86(8)  =  .  999 

2415 

86(9)=. 999999 

2420 

86 ( 1 0 ) = 1 

2425 

B6( 1 )=. 0266 

' . 0372 

2430 

B6(2)=. 0723 

'  .  15 

2435 

B6 ( 3 ) = . 339 

!  .  474 

2440 

B6(4)=. 765 

!  .  9 

2445 

B6 ( 5 ) = 1 . 1933 

00 

2450 

B6(6)=l . 717 

2455 

B6(7)=2. 5005 

!  3. 662 

2460 

B6(8)=4. 0335 

'6.215 

2465 

B6C9)=5. 538 

'8.84 

2470 

B6( 10)=  10 

2475 

RETURN 
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2480  !  *********START  OF  INVERSE  PROCEDURE***** 

2485  '  THIS  PROGRAM  IS  CALLED  INVELLCOV.  <X,Y,Z)  IS  A  POINT  IN  A  CARTESIAN  CO- 

2486  i  ORDINATE  SYSTEM.  INVELLCOV  RETURNS  THE  RADIUS  OF  THE  SPHERE  WITH  CENTER 

2490  !  ( H  x , H  y , H  z  >  WHICH  HAS  P3  OF  THE  NORMAL  ELLIPSOIDAL  DISTRIBUTION  WITH 

2491  !  MEAN  (0,0,0)  AND  STANDARD  DEVIATIONS  Sx,Sy,Sz  IN  THE  X,Y,Z  DIRECTIONS, 

2495  '  RESPECTIVELY.  THE  RADIUS  R  IS  GENERALLY  GIVEN  CORRECTLY  TO  AT  LEAST  6- 

2496  '  DECIMAL  DIGITS.  THE  INPUT  P3  SHOULD  SATISFY  ABS  (  P3- 1 /2  X.  499999  .  IF 

2500  '  P3  IS  IN  [0,1)  BUT  THE  ABOVE  INEQUALITY  IS  NOT  SATISFIED,  THEN  THE 

2501  1  OUTPUT  R  MAY  NOT  BE  CORRECT  TO  6-DECIMAL  DIGITS.  A  MAXIMUM  OF  30  ITERA- 

2505  !  TIONS  IS  ALLOWED.  IF  E5=l  THEN  30  ITERATIONS  HAVE  OCCURRED.  IF  E4=l 

2506  1  THEN  P3  IS  IN  (0,1)  BUT  NOT  IN  C 1 ( -6 ) , . 999999 3 .  IF  P3  IS  NOT  IN  [0,1] 

2510  !  THEN  A  VALUE  OF  -IE-99  IS  RETURNED  FOR  R.  SUBROUTINES  USED  ARE  ELLCOV, 

2511  !  ELLRC ,  AND  GAMINV. 

2515  DEF  FNInvel 1 cov<P3, Hx, Hy , Hz, Sx, Sy, Sz, E4, E5> 

2520  COM  X<43) , Y<43) , A6< 10) , B6< 10) , Q9, V, A7, P4, Ve, T6, PI , T9 
2525  K0=Hx*Hx+Hy*Hy+Hz*Hz 

2530  K2=SGR(K0) 

2535  S=MAX(Sx, Sy, Sz) 

2540  A4=. 797884560803  !  SQR(2/PI) 

2545  R3=E5=E4=0 

2550  W4=Sx*Sx+Sy*Sy+Sz*Sz 

2555  IF  ABS(P3-.  5)0. 499999  THEN  2595 

2560  IF  ABS(P3-. 5)<=. 5  THEN  2575 

2565  R3=- 1 E99  !****"P3  IS  UNACCEPTABLE" 

2570  RETURN  R3 

2575  E4= 1  !****P3  NOT  IN  [ 1 E-6 ,. 999999 ] . 

2580  IF  P3=0  THEN  RETURN  0 

2585  IF  P3< 1  THEN  2595 

2590  RETURN  1E99 

2595  Q9=0 

2600  X8=. 0000005 

2605  IF  P3>. 00001  THEN  2620 

2610  X8  =  5E-  1 0 

2615  GOTO  2650 

2620  IF  P3>.001  THEN  2635 

2625  X8=. 000000025 

2630  GOTO  2650 

2635  IF  P3<. 999995  THEN  2650 

2640  X8=. 00000005 

2645  IF  P3>. 9999975  THEN  X8=lE-8  ! NEW 

2650  P5=P3*S/A4 

2655  Q4=MAX(3*P3*Sx»Sy*SzxA4, P5*P5*P5) 

2660  P4=0 

2665  IF  P3  >  = . 5  THEN  2685 

2670  C2=10*X8*< 1-P3) 

2675  R3  =  Q4Xl/3) 

2680  GOTO  2700 

2685  R3=K2 

2690  IF  Q 4  > K 2 * K 2 * K. 2  THEN  R3  =  G4Xl/3> 

2695  C2= 1 0*X8*P3 

2700  Rmin=Q4=R3 

2705  P4=FNE1 1 cov(R3, Hx, Hy, Hz, Sx, Sy, Sz, K0, K2, S) 

2710  IF  P4  >  =  P3  THEN  RETURN  R3 
2715  P5=P4-P3 

2720  FOR  J=1  TO  10 

2725  IF  P3<  =A6 ( J )  THEN  2735 
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I 

5: 


ft 

ft 


3H 


2730 

2735 

2740 

2745 

2750 

2755 

2760 

2765 

2770 

2775 

2780 

2785 

2790 

2795 

2800 

2805 

2810 

2815 

2820 

2825 

2830 

2835 

2840 

2845 

2850 

2855 

2860 

2865 

2870 

2875 

2880 

2885 

2890 

2895 

2900 

2905 

2910 

2915 

2920 

2925 

2930 

2935 

2940 

2945 

2950 

2955 

2960 

2965 

2970 

2975 

2980 

2985 

2990 


! l+<Hx*Hx+Hy*Hy+Hz*Hz)/W4 


NEXT  J 

R3=K2+B6< J)*S 
Rmax=Q5=R3 
T  6=0 

P4  =  FNE 1 1 cov<R3, Hx, Hy , Hz, Sx, Sy , Sz, K0, K2, S> 

IF  P4<  =P3  THEN  RETURN  R3 
P6=P4-P3 
1 6  =  Q4 
I  7  =  Q5 

!  *************GRUBB'S  ESTIMATE  FOR  R3***************** 

D2  =  0 
D3»l .  1 

IF  P3<  = . 8  THEN  2810 
D3- 1 . 25 

IF  P3<  = . 9  THEN  2810 

D3=  1 . 9 

D4=P4=P3 

M0=  1  +K0/W4  !  l+<Hx*Hx+Hy*Hy+Hz*Hz)/W4 

Wx=Sx*Sx/W4 

Wy=Sy*Sy/W4 

Wz=Sz*Sz/W4 

Gx=Hx/Sx 

Gy=Hy/Sy 

Gz=Hz/Sz 

V4  =  2*  (.  Wx*Wx*  ( 1  +2*Gx*Gx)+Wy*Wy*<l+2*Gy*Gy)+Wz*Wz*(  l+2*Gz*Gz>  ) 

V5  =  8*<Wx*Wx*Wx*<l+3*Gx*Gx)+Wy*Wy*Wy*<:i+3*Gy*Gy>+Wz*Wz*Wz*<:i+3*Gz*Gz>> 

W2=.5*V5/V4 

V5=V5*V5/< V4*V4*V4> 

R4  =  R3 

IF  P4= 1  THEN  P4=. 999999999 

R3  =  W4*<<FNGaminv<4''V5,X,0,P4,  1-P4,  Ierr)-4/V5>*W2+M0> 

IF  <R3<=Q4*Q4)  OR  <R3>=Q5*Q5>  THEN  3090 
R3=SQR ( R3 ) 

IF  ABS<R3-R4X5E-9*R3  THEN  3090 

P4  =  FNE1 1  com  <R3, Hx , Hy , Hz , Sx, Sy , Sz , K0 , K2, S> 

1 8  =  P4 

IF  P4  >0  THEN  Q9  =  Q9+ 1 

IF  P4 >5E-9  THEN  2940 

IF  R3<  =Q4  THEN  3090 

Q4=R3 

P5=P4-P3 

GOTO  3090 

W5=P4-P3 

IF  ABS ( W5  X  =X8  THEN  RETURN  R3 
IF  W5<  0  THEN  3020 
Q5  =  R3 
1 7*Q5 
P6  =  W5 

IF  D2< 0  THEN  3175 

IF  Q9 >2  THEN  3090 

IF  1 8  > 1 0*P3  THEN  3090 

IF  <  P3  > . 0 1 >  OR  <  P  4  = 1 >  THEN  3005 

IF  ABS  <  W5 ) > . 1 *P3  THEN  3005 


,  •  .  .. 


2995 

3000 

3005 

3010 

3015 

3020 

3025 

3030 

3035 

3040 

3045 

3050 

3055 

3060 

3065 

3070 

3075 

3080 

3085 

3090 

3095 

3100 

3105 

3110 

3115 

3120 

3125 

3130 

3135 

3140 

3145 

3150 

3155 

3160 

3165 

3170 

3175 

3180 

3185 

3190 

3195 

3200 

3205 

3210 

3215 

3220 

3225 

3230 

3235 

3240 

3245 

3250 
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P4=D4=P3-2*W5 
GOTO  3010 
P4=D4=D4~D3 
D2=  1 

GOTO  2870 
Q4  =  R3 
1 6  =  Q4 
P5  =  W5 
D5=2-D3 

IF  D2 >0  THEN  3175 

IF  <  Q9  >4 )  OR  < I8< . 01*P3>  THEN  3090 

IF  ( P3 > . 0 1 )  OR  ( P4  =  0 )  THEN  3070 

IF  ABS  <  W5 ) > . 1  *P3  THEN  3070 

P4=D4=P3-2*W5 

GOTO  3075 

P4=D4=D4-D5 

D2  =  -  1 

GOTO  2870 

!  *************ESTIMATE  FOR  R3  USING  CUBE*************** 
PI  1=P1=FNE1 1 rc  C H x , Hy, Hz, Sx, S y , Sz, P3, 0, G  5 , T  6  > 

R4  =  R3 

IF  T 6 <  =  V e  THEN  3110  !  USED  FOR  MAX.  OF  T6. 

Ve=T6  !  USED  FOR  MAX.  OF  T6. 


IF  T 6 <  =  V e  THEN  3110  !  USED  FOR  MAX.  OF  T6. 

Ve=T6  !  USED  FOR  MAX.  OF  T6. 

IF  P  1  <  =  Q  4  THEN  3130 
R3  =  P  1 

GOSUB  3740 

IF  P4>=0  THEN  RETURN  PI 
P 1 = 1 . 73205 1 *P 1 
IF  P 1 <  Q5  THEN  G5  =  P1 

IF  <  P3<  =  . 2 )  AND  <  P 1 1 =  Q  4  >  THEN  3160 
R4  =  R3 

R3=.  5*<Q4  +  Q5'j 

GOSUB  3740  !  COMPUTES  ELLC0V(R3>  AND  MAKES  PROPER  STORAGES. 

IF  ABS<P4X=X8  THEN  3660 

IF  <  1 8<  =5E-9 )  OR  <ABS<Q5-Q4>>. 1*R3>  THEN  3145 
IF  <ABS(R3-R4)<=. 005*R3)  OR  < ABS ( P6-P5 ) < = . 00 1 *P3 >  THEN  3190 
R4  =  R3 

R3=Q4-<Q4-Q5)/CP5-P6)*P5  !  REGULA-FALSI 

GOTO  3200 
R4  =  R3 

R3= . 5*(Q4+G5> 

GOSUB  3740 

IF  ABS < P4 ) <  =X8  THEN  3660 
IF  ABS(Q5-Q4><=. 000001*R3  THEN  3315 

!  ********** ******STEPP I NG  PROCEDURE**************** 

I  9= . 2  *  <  Q5-Q4  > 

P7  =  P5 
R3  =  Q4 

IF  ( P3< . 2 )  OR  (P6+P5>0>  THEN  3255 
I  9  =  -  I  9 
P7  =  P6 
R3  =  Q5 


m 
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FOR  J2= 1  TO  4 
R4=R3 
R3=R3+  1 9 
T9=3270 
GOSUB  3740 

IF  flBS <  P4  )  <  =  X8  THEN  3660 
IF  Q9  >30  THEN  3725 
IF  P4*P7 >0  THEN  3305 

IF  < P3< . 2  )  RND  <RBS< I8/P3-5. 05) >4. 95)  THEN  3220 
GOTO  3310 
NEXT  J2 

IF  <  A  B  S  ( 17-16) >.  005*R3)  OR  < RBS < P6-P5 ) > . 00 1 *P3 )  OR  <I8> 
R3= . 5*  C  Q4  +  Q5 )  !  R3= < 2*Q4+Q5 ) /2 

GOTO  3345 

!  ♦*********#*****KING/S  PROCEDURE****************** 

R4  =  R3 

R3  =  Q5-P6*<Q5-Q4)-'<P6-P5) 

T9  =  3400 

IF  Q9>=30  THEN  3725 

P4  =  FNE 1  1  cov<R3, Hx, Hy, Hz, Sx, Sy, Sz, K0,K2, S) 

IF  P4<  >0  THEN  Q9  =  Q9+ 1 
P4=P4-P3 

IF  RBS  < P4 ) <  =X8  THEN  3660 
IF  P4  >0  THEN  3390 

1 6  =  R3 

I  3  =  -P4 
GOTO  3400 

1 7  =  R3 
I  4  =  P4 

IF  I 7<  =  I  6  THEN  3555 

IF  RBS< I7-I6)<=. 0000005*R3  THEN  3660 

IF  <RBS(R4-R3)<  =  . 000001*R3)  AND  < RBS < P4 ) < =C2 )  THEN  3660 

IF  Q9 >20  THEN  3605 

IF  P4*P6  >0  THEN  3455 

R6  =  Q5 

Q5  =  Q4 

Q4  =  R6 

R6  =  P5 

P5  =  P6 

P6  =  R6 

P5=P5*P6/<P6+P4) 

P6  =  P4 
Q5  =  R3 
R4  =  R3 

R3=Q5-(Q4-Q5)*P6/(P5-P6) 

T9=3480 

IF  Q9>=30  THEN  3725 

P4  =  FNE  1  1 c  ov  <  R3 , Hx, Hy, Hz, Sx, Sy , Sz, K0 , K2 , S  ) 

IF  P4<  >0  THEN  Q9  =  Q9+ 1 
P4=P4-P3 

IF  <RBS<P4X=X8>  RND  <P5<>-P3)  AND  CP601-P3)  THEN  3660 
IF  P4  >0  THEN  3530 
I  6  =  R3 
1 3  =  -P4 
GOTO  3540 
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3645 
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1 7  =  R3 
1 4  =  P4 

IF  ABS< I7-I6><=. 0000005*R3  THEN  3660 
IF  Q9 >20  THEN  3605 
IF  1 7  > 1 6  THEN  3570 
R3=  1 6 

IF  ABS<  I4XABSCI3)  THEN  R3=I7 
RETURN  R3 

IF  <ABS<R4-R3)<=. 000001*R3>  AND  < ABS C P4 > < =C2 >  THEN  3660 

IF  P4*P6 >  =  0  THEN  3455 

Q4  =  Q5 

Q5  =  R3 

P5=P6 

P6  =  P4 

GOTO  3330 

Q4=  I  6 

Q5=  1 7 

R4  =  R3 

R3= . 5*  <  Q4  +  Q5 ) 

T9=3685 
GOSUB  3740 

IF  ABS<Q5-Q4)<=. 000001*R3  THEN  3660 
IF  ABS ( P4 ) <  =X8  THEN  3660 
IF  Q9  >30  THEN  3725 
GOTO  3615 

!  ************ ****EX ITS******** *****##*8 
IF  P4*P6<  0  THEN  3680 
IF  P4*P5  >0  THEN  3700 
R3=R3-<Q4-R3)*P4/<P5-P4> 

RETURN  R3 

R3=R3-<Q5-R3)*P4/<P6-P4) 

RETURN  R3 
R3= . 5* ( Q4+Q5 ) 

RETURN  R3 

IF  P4< 0  THEN  3715 

R3  =  M I N  <  R3 , R4 ) 

RETURN  R3 
R3  =  MAX  <  R3 , R4  > 

RETURN  R3 

E5= 1  !R3  GIVEN  AFTER  30  ITERATIONS. 

RETURN  R3 

!  *********  A  SUBROUTINE  FOR  P4************** 

P4=FNE 1 1 cov<R3, Hx, Hy, Hz, Sx, Sy, Sz,K0, K2, S> 

IF  P4  >0  THEN  Q9=Q9+ 1 

1 8  =  P4 
P4=P4-P3 

IF  P4  >0  THEN  3790 

Q4  =  R3 

1 6  =  R3 

P5  =  P4 

I 3=-P4 

RETURN 
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3790 

3795 

3800 

3805 

3810 

3815 

3820 

3825 

3826 

3830 

3831 

3835 

3836 

3837 

3840 

3841 
3845 
3850 
3855 
3860 
3865 
3870 
3875 
3880 
3885 
3890 
3895 
3900 
3905 
3910 
3915 
3920 
3925 
3930 
3935 
3940 
3945 
3950 
3955 
3960 
3965 
3970 
3975 
3980 
3985 
3990 
3995 
4000 
4005 
4010 
4015 
4020 
4025 
4030 


Q5=R3 
1 7=R3 
P6  =  P4 
1 4  =  P4 
RETURN 
FNEND 

!  *********************ELLCOV******************************* 

THIS  PROGRAM  IS  CALLED  "ELLCOV“.  (X,Y,Z)  IS  AN  ELEMENT  OF  A  CARTESIAN 
COORDINATE  SYSTEM.  ELLCOV  RETURNS  THE  PROBABILITY  OF  A  SHOT  FALLING  UNDER 
AN  ELLIPSOIDAL  NORMAL  DISTRIBUTION,  IN  A  SPHERE  WITH  CENTER  (Hx,Hy,Hz) 

AND  RADIUS  R3.  THE  DISTRIBUTION  HAS  MEAN  2ERO  AND  STANDARD  DEVIATIONS  Sx, 
Sy,Sz  IN  THE  X , Y , Z  DIRECTIONS,  RESPECTIVELY.  THE  INPUT  PARAMETERS  ARE 
R3,Hx,Hy,Hz,Sx,Sy,Sz,K0,K2,S,  WHERE  K0=Hx*Hx+Hy*Hy+Hz*Hz ,  K2=SQR(K0) 

AND  S=MAX(Sx, Sy, Sz ) . 

THE  OUTPUT  P  IS  GENERALLY  ACCURATE  TO  AT  LEAST  6-DECIMAL  DIGITS. 
SUBROUTINES  USED  ARE  CIRCV,  DAWS,  DXDAWS,  DXERF ,  DXERF3 ,  ELLCV ,  EQSIG, 
ERF,  ERF3 ,  ERFC,  ERFC1 ,  FN. 

DEF  FNE  1  lco<vKR3,Hx,Hy,Hz,Sx,Sy,Sz,K0,K2,S> 

COM  X(43) , Y(43) , A6( 10) , B6( 10) , Q9, V, A7,P4, Ve, T6, PI , T9 

C5=6. 02738  15.6123 

C6=. 398942280402  !1/SQR(2*PI) 

C7=. 707106781 187  ! SQR  < . 5 ) 

C8=l. 41421356237 
F2=4 . 262 

C10=. 999999998334  ! ERF  <  F2  ) 

Bl=. 564189583548  !1/SQR(P1) 

A4=. 797884560803  !SQR<2/PI) 

Z4=. 000000025 
H  =  ABS  <  Hx  ) 

K=ABS ( Hy ) 

L  =  ABS  <  Hz ) 

S 1  =Sx 
S2=Sy 
S3  =  Sz 
P4  =  0 
T5  =  0 

Hl=1.5*Sx*Sy*Sz*Z4xC6 

K7=R3*R3*R3 

IF  K7<  =H 1  THEN  RETURN  P4 

A3=<K2-R3)*C7/S 

A5=EXP(-R3*A3) 

IF  K2< R3  THEN  3985 
IF  K7*A5<  =H 1  THEN  RETURN  P4 
GOTO  4010 

P4=FNEqsi g<R3, K0, K2, S, A3, A5, A4, C7>  !  ELLCOV  FOR  EQUAL  SIGMAS 

V  =  -3 

IF  P4< 1 -Z4  THEN  4010 
P4=  1 

RETURN  P4 

IF  (SxOSy)  AND  (SxOSz)  AND  (SyOSz)  THEN  4305 
IF  (SxOSy)  OR  (SyOSz)  THEN  4035 
IF  R3  >K2  THEN  RETURN  P4 
P4=FNEqsig(R3,K0,K2,S,A3,A5,A4,C7) 

RETURN  P4 
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IF  Sx  =  Sy  THEN  4090 
IF  Sx=Sz  THEN  4070 
S3  =  Sx 

5 1  =Sz 
L=RBS ( Hx ) 

H=RBS ( Hz ) 

GOTO  4090 

52  =  Sz 

53  =  Sy 

K  =  FIBSc:Hz) 

L  =  F)BS  <  Hy  ) 

IF  H  +  KO0  THEN  4280 

!  ************31=82  RND  H=K=0********************** 

Ll=flBS<  <Sl*Sl-S3*S3)x<Sl*Sl  )  > 

L2=SQRCL1 ) 

L4=<L+R3)*C7xS3 
L5«<L-R3)*C7/'S3 
IF  L5  >F2  THEN  RETURN  P4 
W9  =  -  1 

J3=EXP<-L5*L5> 

S0=C7/<S3*L2) 

IF  L=0  THEN  4240 
IF  R3< L  THEN  4160 
P4=FNErf <L4)-FNErf3(L5, J3) 

GOTO  4170 

W9=EXP<-2*R3*L/<S3*S3>) 

P4=J3*<FNErfc 1 0.5, 1 ) -W9*FNErf c 1 < L4 , 1)) 

F1=<L+L1*R3)*S0 
S0=CL-L1-*R3)*S0 
IF  S 1 <  S3  THEN  4220 

V  =  -5 

IF  S  0  <  0  THEN  4210 

IF  W9<  0  THEN  UI9  =  EXP<-2*L*R3/<S3*S3>  ) 

P4=. 5*<P4-J3/L2*<FNErfc 1 < S0 , 1 > -W9*FNEr f c 1 < F 1 , 1 >  >  ) 

RETURN  P4 

P4=. 5*<P4-l/L2*EXPO. 5*<R3*R3-L*L/L1 >/<Sl*Sl >  >*<FNErf <F1 >-FNErf <S0>  )  ) 
RETURN  P4 

V  =  -6 

IF  W9<  0  THEN  W9  =  EXP(-2*L*R3/(S3*S3) > 

P4=.  5*P4-B1 * J3/L2* < FNDaws (FI ) -W9*FNDaws (  S0  ) > 

RETURN  P4 

IF  S3 >S 1  THEN  4260 

V  =  -5 . 5 

P4=FNErf3(L4, J3)-L4*EXP<-. 5*R3*R3x < S 1 *S 1 ) ) *FNDxerf 3 < L4*L2 , - 1 > 

RETURN  P4 

V  =  -6 . 5 

P4=FNErf3<L4, J3 > -2*B 1  * J3*L4*FNDxdaws ( L4*L2 > 

RETURN  P4 

!  ******************Sl=S2***** ************* 

V=  1 

D  =  SQR  <  H*H  +  K*K> 

GOTO  4415 
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!  *****************H=K=0******************** 

IF  ( H  +  K  =  0 )  AND  (ABSCS1/S2-50000. 000005K49999. 999995)  THEN  4385 
IF  <H  +  L  =  0)  AND  <ABS(Sl/$3-50000. 000005X49999. 999995)  THEN  4355 
IF  (K+LO0)  OR  (ABS(S2/S3-50000. 000005)  >49999.  999995)  THEN  4510 
L 1  =S  1 

5 1  =S3 
S3  =  L  1 
L 1  =H 
H  =  L 
L=L  1 

GOTO  4385 
L 1  =S2 

52  =  S3 

53  =  L  1 
L 1  =  K 
K  =  L 

L  =  L  1 

IF  S 1 >S2  THEN  4405 
W1=S1 

5 1  =S2 

52  =  W  1 
D=S2/S1 
V  =  0 

X1=C7/S3 
A5=X  1  *  (  L-R3  ) 

IF  A5 >F2  THEN  RETURN  P4 
W9=  1 

IF  A5  >-F2  THEN  4475 
A5=-F2 
B5  =  F2 
X2=L*X 1 

IF  X2 >  =  F2  THEN  4825 

W9=0 

B5  =  X2 

GOTO  4825 

B5=X 1  *L 

IF  B5  >F2  THEN  4495 
W9  =  0 

GOTO  4825 
B5  =  F2 
GOTO  4825 

!  ***************GENERAI_  CASE************************ 

V=2 

HC 1 )=H 
H  <  2 ) =K 
H  <  3 ) =L 
S< 1 )=S1 
S ( 2 ) =S2 
S ' 3 ) =S3 

!  **************F IX  ORDER  INTEGRATION*************** 

T (  1 )  «MAX  <  S 1 , S2 ) 

T<2)=MAX<S2, S3) 

T  <  3 ) *MAX  <  S3 , S 1 ) 
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4565  W8U)=S1/S2 

4570  IF  W8  < 1 >  > 1  THEN  W8 < 1 )  =  1 /W8 < 1  ) 

4575  WS<2)=S2/S3 

4580  IF  W8  <  2 ) > 1  THEN  W8 ( 2 > = 1 /W8 < 2 > 

4585  W8<3)=S3/S1 

4590  IF  W8 < 3 ) > 1  THEN  W8 < 3 ; -  1 /W8 < 3 > 

4595  J=3 

4600  IF  W8(1>>W8<2)  THEN  4645 

4605  IF  W8(1)<>W8(2)  THEN  4620 

4610  IF  T(1XT(2)  THEN  J=1 

4615  GOTO  4670 

4620  J=1 

4625  IF  Ul 8  <  2  >  >W8  ( 3 )  THEN  4670 
4630  IF  W8  <  2  X  >W8 < 3 )  THEN  4660 
4635  IF  T<2XT<3>  THEN  J  =  2 
4640  GOTO  4670 

4645  IF  M8 < 1 >  >W8 < 3 )  THEN  4670 
4650  IF  W8<1X>W8<3>  THEN  4660 

4655  IF  T<1)>T<3)  THEN  4670 

4660  J  =  2 

4665  !  ***********FIND  LIMITS  OF  INTEGR8TI0N************** 

4670  X1=C7/S<J) 

4675  85< J)=<H< J)-R3)*X1 

4680  IF  85  ( J  ) >=F2  THEN  RETURN  P4 

4685  W8  <  J )  =  1 

4690  IF  85  <  J ) >-F2  THEN  4730 

4695  85  <  J ) =-F2 

4700  B5  <  J ) =F2 

4705  X2  =  H  <  J ) *X 1 

4710  IF  X2  >  =  F2  THEN  4755 

4715  B5  <  J  ) =X2 

4720  W8 ( J ) =0 

4725  GOTO  4755 

4730  B5< J)=H< J)*X1 

4735  IF  B5 <  J ) >F2  THEN  4750 

4740  W8  <  J ) =0 

4745  GOTO  4755 

4750  B5(J)=F2 

4755  IF  J=3  THEN  4810 

4760  IF  J=2  THEN  4790 

4765  H  =  H  <  3 ) 

4770  L  =  H  <  l  > 

4775  S 1 =S  <  3  > 

4780  S3  =  S  (  1 ) 

4785  GOTO  4810 

4790  L  =  H (2) 

4795  S3=S ( 2 > 

4800  K  =  H  <  3 ) 

4805  S2  =  S  <  3 ) 

4810  85  =  85  <  J  > 

4815  B5  =  B5  <  J  ) 

4820  W9=W8<J> 
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!  ********** **SPEC I fll_  CASE************* 

Py  =  0 

IF  A5 >-F2  THEN  4915 

K7=L-C8*S3*B5 

T8=L-C8*S3*A5 

K7=R3*R3-K7*K7 

T8=R3*R3-T8*T8 

IF  <K7<0>  OR  < T8< 0 )  THEN  4915 
K7=SQR  <K7  ) 

T8=SQR<T8> 

IF  ABSCT8-K7) >K7*lE-8  THEN  4915 
V  =  -9  +  V 

IF  VO-7  THEN  4895 

P4  =  C 1 0*FNE 1 lcv<. 5*<K7+T8) , H, K, SI ,  S2> 

RETURN  P4 

P4  =  C10*FNCi rcv< . 5* ( K7+T8 > , D ,  S 1 , V ) 

RETURN  P4 

!  ******** *PT ART  OF  GAUSSIAN  INTEGRATION*********** 

K7  =  MAX  <  S 1 ,  S2 ) 

T8=C5*K7+SGR<H*H+K*K> 

T8=R3*R3-T8*T8 
IF  T8<  0  THEN  4985 
A3  =  SQR  C  T8 ) 

T8=<L-A3>/(C8*S3) 

IF  T  8  >  =  B5  THEN  4985 
B5=T8 

K7=<L+A3)/(C8*S3) 

IF  K7<  4 . 4  THEN  4975 
Py  =  P4=. 5*FNErfc  <T8) 

GOTO  4985 

Py=P4=. 5*<FNErf <K7)-FNErf <T8) ) 

IF  <ABS<B5-A5X=ABS<B5  +  A5)*1E-10>  OR  <B5<=A5>  THEN  RETURN  P4 
A3  =  A5 

FOR  I =- 1 2  TO  12 
IF  1=0  THEN  5020 

T5=. 5*( <B5-A5)*SGN< I ) *X < 3 1 +ABS < I ) ) +B5+ A5 > 

X2  =  FNF/i<T5,  R3,  H,  K,  L,  SI , S2, S3, V, D, W9> 

IF  X2 > IE- 1 0  THEN  5025 
A3  =  T5 
NEXT  I 

IF  X2< 1 E-7  THEN  5115 
A5=A8= A3 

IF  <T5-A5)*X2<=lE-9  THEN  5115 
T8  =  A5 

K7=<T5-A5>*< 1-XC43) )  ! . 5*K7= . 5* < T5+A5- ( T5-A5 > *X C 43 ) ) -A5 

IF  <ABS<K7X=ABSCA5>*1E-10>  OR  <K7<=0>  THEN  5115 

T8=T8+K7 

X3  =  FNFr,<T8,  R3,  H,  K,  L,  SI  ,  S2,  S3,  V,  D,  W9>  !  COMPUTATION  OF  INTEGRAND 

IF  X3 >1E- 1 0  THEN  5080 

A8  =  T8 

GOTO  5055 

A3=A5=A8 

IF  X3<  *  IE-7  THEN  5115 
IF  <T5-A5)*X2<=lE-9  THEN  5115 


m 
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5095  !  ************GAUSSI AN  RULE.  Pz  GIVES  THE  I NTEGR AL  FROM 

5096  !  A5  TO  T5  OF  THE  INTEGRAND  F  (X>=FNFn<:x> .  *********** 

5100  CALL  Gquad ( A5 , T5 , I  err, 1 , Pz , E , 1 2 , N2 , R3 , H , K , L , S 1 , S2 , S3 , V, D , W9 > 

5105  A5=T5 

5110  GOTO  5130 

5115  Pz  =  0 

5120  A5  =  A3 

5125  !  ************GAUSSIAN  RULE.  P4  GIVES  THE  INTEGRAL  FROM 

5126  !  A5  TO  B5  OF  THE  INTEGRAND  F(X)=FNFn(X) . *********** 

5130  CALL  Gquad < A5, B5, Ierr ,  1 ,  P4,  E,  1 2,  N2,  R3,  H,  K,  L,  SI ,  S2,  S3,  V,  D,  W9  > 

5135  P4=B1*<P4+Pz)+Py 

5140  IF  P4  > 1  THEN  P4=l 

5145  RETURN  P4 

5150  FNEND 

5155  !  ******************EQSIG************************************** 

5160  !  FNEqsig  COMPUTES  ELLCOV  FOR  EQUAL  SIGMAS. 

5165  !  SUBROUTINES  NEEDED:  ERF,  ERF3,  ERFC1,  REXP. 

5170  DEF  FNEqsi g<R3, K0, K2, S, A3, A5, A4, C7> 

5175  IF  K0O0  THEN  5195 

5180  V=-3 

5185  P4  =  FNErf3<-A3,  A5>-A4*R3*A5/'S 

5190  RETURN  P4 

5195  K1=(R3+K2>*C7/S 

5200  P5=2*R3*K2/(S*S) 

5205  Q5=FNRexp<-P5) 

5210  IF  R3 >K2  THEN  5230 

5215  Q4  =  FNErfc  1  <K1  ,  1  > 

5220  P4  =  A5*  <  FNEr f  c 1 <A3, 1 ) -Q4+P5*Q5*Q4 > 

5225  GOTO  5235 

5230  P4=FNErf <K1 >-FNErf3<A3, A5) 

5235  P4=. 5*P4-A4*A5*R3*Q5/S 

5240  V  =  -3 . 5 

5245  RETURN  P4 

5250  FNEND 

5255  !  *******  ***************GAM  I  NV***********************  *********** 

5260  !  SUBPROGRAM  FUNCTION  GAM  I N V < A , X , X0 , P , Q ,  I  ERR ) 

5265  !  INVERSE  INCOMPLETE  GAMMA  RATIO  FUNCTION 

5270  !  GIVEN  A>0  AND  NONNEGATIVE  P  AND  Q  WHERE  P+Q=l,  THEN  X  IS  OBTAINED 

5271  !  WHERE  P(A,X)=P  AND  Q<A,X)=Q. 

5275  !  SCHRODER  ITERATION  IS  EMPLOYED.  THE  ROUTINE  ATTEMPTS  TO  COMPUTE  X 

5276  !  TO  10  SIGNIFICANT  DIGITS. 

5280  !  X  IS  A  VARIABLE.  IF  P=0  THEN  X  IS  ASSIGNED  THE  VALUE  0,  AND  IF  Q=0 

5281  !  THEN  X  IS  SET  TO  THE  LARGEST  FLOATING  POINT  NUMBER  AVAILABLE. 

5285  !  OTHERWISE,  GAM  I NV  ATTEMPTS  TO  OBTAIN  A  SOLUTION  FOR  P(A,X)=P  AND 

5286  !  GKR,X)=Q.  IF  THE  ROUTINE  IS  SUCCESSFUL  THE  SOLUTION  IS  IN  X. 

5290  !  X0  IS  AN  OPTIONAL  INITIAL  APPROXIMATION  FOR  X.  IF  THE  USER  DOES  NOT 

5291  !  WISH  TO  GIVE  AN  INITIAL  APPROXIMATION,  THEN  SET  X0  <=  0. 

5295  !  Ierr  IS  A  VARIABLE  THAT  REPORTS  THE  STATUS  OF  THE  RESULTS.  WHEN  THE 

5296  !  ROUTINE  TERMINATES,  Ierr  HAS  ONE  OF  THE  FOLLOWING  VALUES: 

5300  !  I  err  =  0  THE  SOLUTION  WAS  OBTAINED.  ITERATION  WAS  NOT  USED. 

5305  !  I er r  >0  THE  SOLUTION  WAS  OBTAINED.  Ierr  ITERATIONS  WERE  DONE. 
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3310 

1 

I err  =  -2  (INPUT  ERROR)  Fl<=0. 

5315 

! 

I err=-3  X  IS  TOO  SMALL.  X  =  0  GIVEN. 

3320 

! 

I err=-4  (INPUT  ERROR)  P+Q  #  1. 

3325 

i 

I err=-6  20  ITERATIONS  DONE.  MOST  RECENT  VALUE  OF  X  GIVEN. 

5326 

! 

CANNOT  OCCUR  IF  X0  <=  0. 

5330 

j 

I e r r  =  - 7  ITERATION  FAILED.  NO  VALUE  FOR  X  IS  GIVEN.  THIS  MAY 

5331 

( 

OCCUR  WHEN  X  IS  APPROXIMATELY  0. 

5335 

! 

I e r r  =  - 8  A  VALUE  FOR  X  HAS  SEEN  OBTAINED,  BUT  GAM I N V  I 

S  NOT 

5336 

! 

CERTAIN  OF  ITS  ACCURACY.  ITERATION  CANNOT  BE 

DONE 

5340 

! 

IN  THIS  CASE.  IF  X0  <=  0,  THIS  CAN  OCCUR  ONLY 

WHEN  P 

5341 

j 

OR  Q  IS  APPROXIMATELY  0.  IF  X0  >  0  THEN  THIS 

CAN  OCCUR 

5345 

i 

WHEN  A  IS  EXCEEDINGLY  CLOSE  TO  X  AND  A  IS  VERY  LARGE 

3346 

! 

(SAY  A  >  1E20). 

5350 

j 

CODE  TAKEN  FROM  THE  FORTRAN  VERSION  BY  AL  MORRIS. 

5355 

i 

SUBROUTINES  NEEDED:  GAMMA I ,  GRATIO,  ALNREL ,  GAMLN,  RCOMP, 

GAMLN 1 

5360 

DEF  FNGaminy(A,X,X0,P,Q,Ierr) 

5365 

C=. 577215664902 

5370 

Ln 1 0=2 . 302585 

5375 

A0=3. 31 125922109 

5380 

A  1  =  1 1 . 6616720289 

5385 

A2=4. 28342155967 

5390 

A3=. 213623493716 

5395 

Bl=6. 61053765625 

5400 

B2=6. 4069159776 

3405 

B3= 1 . 27364489782 

3410 

B4=3. 61 170810188E-2 

5415 

Eps0( 1 )=1E-10 

3420 

Eps0 ( 2 ) = 1 E-8 

5425 

Ami r( 1 )=500 

3430 

Ami n(2)=100 

3435 

Bmi n( 1 )=lE-28 

5440 

Bmin(2)=lE-13 

5445 

Dmi n( 1 )=lE-6 

5450 

Dmin(2)=lE-4 

3455 

Em i n ( 1 ) =2E-3 

5460 

Em i n ( 2 ) =6E - 3 

5465 

To  1 = 1 E-5 

5470 

E= 1 E- 1 2 

3475 

Xm i n= 1 E-99 

5480 

Xmax=9 . 999999999E99 

3485 

X  =  0 

5490 

IF  A<  =0  THEN  6405 

5495 

T=- . 5+P+Q- . 5 

5500 

E2=E+E 

5505 

IF  ABS ( T ) > 1 0*E  THEN  6425 

5310 

I  err  =  0 

5515 

IF  P=0  THEN  RETURN  X 

5320 

IF  Q=0  THEN  6365 

5525 

IF  A=1  THEN  6375 

5530 

Amax=4E- 1 1 /(E*E) 

3535 

I  op=  1 

5540 

IF  E  > 1 E - 1 0  THEN  Iop  =  2 

5545 

Eps=Eps0 ( lop) 

3530 

Xn  =  X0 

5355 

IF  X0>0  THEN  6070 
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5560  !  ****INITAIL  APPROXIMATION  Xn  OF  X  FOR  A<  1  .  **** 

5565  IF  A  > 1  THEN  5795 

5570  G= 1 /FNGammai < . 5+A+ . 5 ) 

5575  Qg  =  Q*G 

5580  IF  Qg=8  THEN  6470 

5585  B  =  Qg/A 

5590  IF  Qg > .  6*R  THEN  5735 

5595  IF  ( A >= . 3 )  OR  <B<.35>  THEN  5620 

5600  T=EXP ( - ( B+C ) > 

5605  U  =  T  *EXP ( T ) 

5610  Xn  =  T  *EXP ( U ) 

5615  GOTO  6070 

5620  IF  B  >= . 45  THEN  5735 

5625  IF  B=0  THEN  6470 

5630  Y=-LOG<B) 

5635  S= . 5- A+ . 5 

5640  Z=L0GCY> 

56^5  T=Y-S*Z 

5650  IF  B< . 15  THEN  5665 

5655  Xn=Y-S*LOG<T)-LOG< l+S/CT+l >  > 

5660  GOTO  6220 

5665  IF  B<  =  . 0 1  THEN  5685 

5670  U=( (T+2*<3-A) )*T+<2-A)*<3-A) >  /  < (T+C5-R) >  *T  +  2  > 

5675  Xn=Y-S*L0G(T>-L0G<U> 

5680  GOTO  6220 

5685  Cl=-S*2 

5690  C2= -S*  < 1 +C 1 ) 

5695  C3  =  S*<  < . 5*Cl+<2-A) )*Cl  +  <2. 5-1 . 5*A> ) 

5700  C4  =  -S*  ( <  <  C 1 '  3+(2. 5-1 . 5* A ) ) *C 1  +  <  <A-6>*R+7> 1 *C  1  +  ( < 1 1  * A-46 > *A  +  47 ) /S > 

5705  C5  =  -S*< < ( C-C1/4-K  1 1 *A- 1 7 ) /6 > *C 1  +  < ( -3*A  + 1 3 > *A- 1 3 > ) *C 1  + . 5* <  <  < 2*A-25 > * A  +  72 > 

*  A  -  6  1 )  )  *  C  1  +  (  (  (25*A-195:,*A  +  477}*A-379)rl2) 

5710  Xn=( < (C5/Y+C4>/Y+C3)/Y+C2)/Y+C1+Y 

5715  IF  A > 1  THEN  6190 

5720  IF  B  >Bfn  i  n  <  I  op  >  THEN  6220 

5725  X  =  Xn 

5730  RETURN  X 

5735  IF  B*Q>lE-8  THEN  5750 

5740  Xn  =  EXP<-(Q/A  +  C)  ) 

5745  GOTO  5770 

5750  IF  P<  = . 9  THEN  5765 

5755  Xn  =  EXP< <FN A  1  nr el < -Q ) +FNGam 1 n 1 < A > > /A > 

5760  GOTO  5770 

5765  Xn=EXP(LOG<P*G>/A) 

5770  IF  Xn=0  THEN  6415 

5775  T=. 5  +  < . 5-Xn/< . 5+A+. 5) ) 

5780  Xn=Xn/T 

5785  GOTO  6070 

5790  i  *****  INITIAL  APPROXIMATION  FOR  A  >  1.***** 

5795  IF  Q<  = . 5  THEN  5810 

5800  W=LOG ( P ) 

5805  GOTO  5815 

5810  W  =  L0G  <  Q ) 

5815  T-SQR ( -2*W ) 

5820  S=T-((<A3*T+A2)*T+A1)*T+A0)/<(<<B4*T+B3)*T+B2)*T+B1)*T+1> 
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5825  IF  Q  > . 5  THEN  S  =  -S 
5830  Rt  a=SQR  C  fl ) 

5835  S2=S*S 

5840  Xn=fi+S*Rta+<S2-l  )  /  3  +  S*  C  S2  -  7  >  ■■■  c  36  *  R  t  a  >  —  (  C3*S2  +  7)*S2-16)-/(  818*fi>+S*C  C  9  *  S  2  + 

256>*S2-433>/(38880*ft*Rta> 

5845  Xn  =  f1fiX<Xn,  0> 

5850  IF  fl<Rmi  n<  lop)  THEN  5870 

5855  X=Xn 

5860  D=.5-Xxfl+.5 

5865  IF  RBS<DX  =  Dmi  n(Iop)  THEN  RETURN  X 
5870  IF  P<  =  . 5  THEN  5940 
5875  IF  X n <  3 * fl  THEN  6220 
5880  Y=-CW+FNGaml  r.  ( fl  >  > 

5885  IF  A<  2  THEN  5920 

5890  GOTO  6220 

5895  D  =  MAX(2, A*CR-1  > > 

5900  IF  Y<Lnl0*B  THEN  5920 

5905  S= . 5-8+ . 5 

5910  Z  =  LOG  <  Y  > 

5915  GOTO  5685 

5920  T=- . 5+R- . 5 

5925  Xn=Y+T*LOG<Xn )-FNfll nr e 1 < - T/ < Xn  + 1 > ) 

5930  X  n  =  Y  +  T  *  L  C  G  X  n  >  -  F  N  fi  1  n  r  e  1  C-T/CXn+l  >  ) 

5935  GOTO  6220 

5940  flpl=fl+l 

5945  IF  Xri>.7*Apl  THEN  6075 

5950  W  =  W  +  F  N  G  am  1  n  (  fi  p  1 

5955  IF  Xri>.15*Apl  THEN  6010 

5960  Rp2=fi+2 

5965  Rp3=A+3 

5970  X  =  EXP((14  +  X)/fl) 

5975  X  =  E  X  P  (  (W  +  X-  LOGC  i+XxRpl*(  1+X/Rp2)  )  )  /  R  > 

5980  X=EXP< (W+X-LOGC 1+X/Rpl*< 1+X/Rp2> ) XR) 

5985  X  =  EXP(  (W  +  X-LOGC  1  +X/flp  1  *  (  1  +X^flp2*  <  1  +X''Ap3  )  )  >  )/R> 

5990  Xn=X 

5995  IF  Xn  > . 0 1 *Ap 1  THEN  6010 

6000  IF  Xn<=Emi n<Iop>*Apl  THEN  RETURN  X 

6005  GOTO  6075 

6010  flpn=Ppl 

6015  T  =  X  n  / '  R  p  n 

6020  S  u  m  =  1  +  T 

6025  flpn=flpn+l 

6030  T=T*Xn/Rpn 

6035  Sum=Sum+T 

6040  IF  T  > 1 E-4  THEN  6025 

6045  T=U-L0G(Sum> 

6050  Xn=EXPC <Xn  +  TXfl) 

6055  Xn=Xn*( 1  - < R*LOG < Xn > -Xn-T > ✓ < fl-Xn > ) 

6060  GOTO  6075 

606 5  *  ****SCHR0DER  ITERATION  USING  P--REF.  V0L.1  OF  HENRICI,  1974,  CP. 529), 

6070  IF  P  > . 5  THEN  6220 

6075  IF  POlE10*X«nin  THEN  6455 

6080  fiml =-. 5+A-. 5 

6085  IF  A<  = Amax  THEN  6100 
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D=.  5  -  X  n  •'  A  +  .  5 
IF  ABS ( D ) <  =E2  THEN  6455 
IF  I  err  >  =  20  THEN  6435 
I err= I  err  + 1 

CALL  Grat i o( A, Xn, Pn, Qn, 0) 

IF  <  Pn  =  0 )  OR  <  Qn  =  0 )  THEN  6455 
R=FNRcomp< A, Xn) 

IF  R=0  THEN  6455 

T=<Pn-P)/R 

W=. 5*(Aml-Xn) 

IF  <ABS<TX  =  .l)  AND  < ABS < W*T ) < = . 1 )  THEN  6165 
X=Xn*< 1-T) 

IF  X<  =0  THEN  6445 
D= ABS ( T ) 

GOTO  6190 
H  =  T *  <  1 +  W*T ) 

X=Xn*< 1-H) 

IF  X<=0  THEN  6445 

IF  <  A  B  S  (  W  ) >=1 >  AND  (ABS(W)*T*T<=Eps>  THEN  RETURN  X 
D=ABS (H) 

Xn  =  X 

IF  D  >To 1  THEN  6085 
IF  D<  =  Eps  THEN  RETURN  X 
IF  ABS<P-Pn)<=Tol*P  THEN  RETURN  X 
GOTO  6085 

!  ****SCHRODER  ITERATION  USING  Q****** 

IF  Q<  =  lE10*Xmi  ri  THEN  6455 

Am  1 =- . 5  +  A- . 5 

IF  A<  =  Amax  THEN  6245 

D  = . 5-Xn/A+ . 5 

IF  ABS ( D ) <  =E2  THEN  6455 

IF  I  err  >  =  20  THEN  6435 

I err= I err+ 1 

CALL  Gr at i o < A , Xn , Pn , Qn , 0 ) 

IF  < Pn  =  0  >  OR  <  Qn  =  0 >  THEN  6455 
R  =  FNRc  omp ( A , Xn ) 

IF  R=0  THEN  6455 
T  =  <  Q-Qn  )  /R 
W=. 5*<Aml-Xn) 

IF  <  ABS  <  T ) <  =  .  1)  AND  < ABS < W*T X  = . 1 )  THEN  6310 
X  =  Xn*  < 1-T) 

IF  X <  =  0  THEN  6445 
D  =  ABS  <  T ) 

GOTO  6335 
H  =  T  *  ( 1 +  W*T ) 

X=Xn*< 1-H) 

IF  X<=0  THEN  6445 

IF  ( ABS  <  W  >  >= 1 )  AND  < ABS < W ) *T*T< =Eps >  THEN  RETURN  X 
D=ABS(H) 

Xn  =  X 

IF  D >To 1  THEN  6230 
IF  D<  =Eps  THEN  RETURN  X 
IF  ABS<Q-QnX=Tol*Q  THEN  RETURN  X 
GOTO  6230 
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6360  !  *****SPEC I AL  CASES******* 

6365  X=Xmax 

6370  RETURN  X 

6375  IF  Q< . 9  THEN  6390 

6380  X=-FNA1 nre  ?  <-P) 

6385  RETURN  X 

6390  X  =  -LOG  <  Q ) 

6395  RETURN  X 

6400  !  *****ERROR  RETURN******** 

6405  I  err  =  -2 

6410  RETURN  X 

6415  I er r=-3 

6420  RETURN  X 

6425  I  err  =  -4 

6430  RETURN  X 

6435  I  err  =  -6 

6440  RETURN  X 

6445  I err--7 

6450  RETURN  X 

6455  X=Xn 

6460  I  err  =  -8 

6465  RETURN  X 

6470  X=Xmax 

6475  I  err =-8 

6480  FNEND 

6485  !  *******************GR AT  10****** ******************** ******** 

6490  REM  SUBROUTINES  NEEDED:  ERF , ERFC1 ,  GAM  1 ,  GAMMA  I ,  RLOG ,  REXP. 
6495  SUB  Grat i o<A, X, Ans, Qans, Ind) 

6500  DIM  A(3),B<3),E<3>,X(3) 

6505  DIM  W<20),D(13),D1<12),D2<10),D3<8>,D4<6>,D5C4),D6(2) 

65.0  E9  =  5E-  1 2 

6515  E  =  5E-  1 3 

6520  A  < 1 )  =  5E-  1 3 

6525  A  <  2 ) =5E-7 

6530  A ( 3 ) =5E-4 

6535  B ( 1 ) =20 

6540  B  <  2  )  =  1 4 

6545  B  <  3 )  =  1 0 

6550  E ( 1 ) =2 . 5E-4 

6555  E  <  2  >  =  . 025 

6560  E  <  3 )  = .  14 

6565  X ( 1 ) =3 1 

6570  X  <  2 )  =  1 7 

6575  X ( 3 ) =9 . 7 

6580  Lnl0=2. 30258509299 

6585  Rpi =. 398942280401 

6590  Spi =1 . 77245385091 

6595  Th=. 333333333333 

6600  D( 1 )=8. 33333333333E-2 

6605  D<2)=-1 . 48148148148E-2 

6610  D <  3  >  =  1 . 15740740741E-3 

6615  D<4)=3. 52733686067E-4 

6620  D<5)=-1 . 78755144033E-4 

6625  D<6>=3. 9 1 9263 1 7852E-5 
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6630  D<7>=-2. 18544851068E-6 

6635  D  <  8  >  =- 1 . 85406221872E-6 

6640  D<9>=8. 29671 134095E-7 

6645  D< 10>=-1 . 76659527368E-7 

6650  D<1 1 >=6. 70785354340E-9 

6655  D< 12)=1 . 02618097842E-8 

6660  D< 13>«-4. 38203601845E-9 

6665  D10=-l . 85185185185E-3 

6670  D1 ( 1 )=-3. 47222222222E-3 

6675  D1 <2>=2. 64550264550E-3 

6680  Dl(3)=-9. 90226337449E-4 

6685  D1 <4)=2. 0576 1 3 1 6872E-4 

6690  D1 C5)=-4. 01877572016E-7 

6695  D1 <6)=-l . 80985503345E-5 

6700  D1 C7)=7. 64916091608E-6 

6705  D1 <8>=-l . 61209008946E-6 

6710  D1 <9>=4. 64712780281E-9 

6715  D1 < 1 0  >  =  1 . 37863344692E-7 

6720  HI ( 1 1 >=-5. 75254560352E-8 

6725  Dl(12)=l. 19516285998E-8 

6730  D20=4. 1 335978836E-3 

6735  D2C 1 )=-2. 68132716049E-3 

6740  D2<2)=7. 71604938272E-4 

6745  D2<3)=2.00938786008E-6 

6750  D2<4)=-1.07366532264E-4 

6755  D2<5)-5.29234488291E-5 

6760  D2(6)=-1.27606351886E-5 

6765  D2<7)=3.4235787341E-8 

6770  D2<8)=1 . 37219573091E-6 

6775  D2<9>=-6. 298992 1 3838E-7 

6780  D2< 10)=1 . 42806142061E-7 

6785  D30=6. 49434 1 56379E-4 

6790  D3< 1 )=2. 29472093621E-4 

6795  D3<2>=-4. 69 1 89494395E-4 

6800  L3<3>=2. 67720632063E-4 

6805  D3(4)=-7. 56180167188E-5 

6810  D3<5)=-2. 3965051 1387E-7 

6815  D3<6)=1 . 10826541 153E-5 

6820  D3<7)=-5. 67495282699E-6 

6825  D3<8)=1 . 42309007324E-6 

6830  D40=-8. 61888290917E-4 

6835  D4< 1 >=7. 8403922 1 72E-4 

6840  D4(2)=-2. 99072480303E-4 

6845  D4<3)=-1 . 46384525788E-6 

6850  D4<4>=6. 64 1 4982 1 547E-5 

6855  D4(5)=-3. 968365047 l 8E-5 

6860  D4(6>=1. 1 3757269707E-5 

6865  D50=-3. 36798553366E-4 

6870  D5< 1 )=-6. 97281375837E-5 

6875  D5<2)=2. 77275324496E-4 

6880  D5<3)=-1 . 99325705162E-4 

6885  D5<4)=6. 79778047794E-5 

6890  D60*5. 3 1 307936464E-4 

6895  L6< 1 )=-5. 92 1 66437354E-4 

6900  D6<2>=2. 70878209672E-4 

6905  D70=3. 44367606892E-4 
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IF  <  A<  0 )  OR  <X<0>  THEN  8055 
IF  <A  =  0>  AND  <X  =  0>  THEN  8055 
IF  A*X=0  THEN  8020 
I=Ind+l 

IF  (lOl)  AND  (102)  THEN  1=3 
Ac  c  =  MAX ( A  < I ) , E) 

E0=E(I> 

X0=X< I > 

*******SELECT  APPROPRIATE  ALGORITHM****** 

IF  A >= 1  THEN  6995 

IF  A= . 5  THEN  7980 

IF  X< 1 . 1  THEN  7370 

T 1 =A*LOG  <  X ) -X 

U=A*EXP ( T 1 ) 

IF  U=0  THEN  8040 
R=U*< 1+FNGaml <A> ) 

GOTO  7605 

IF  A  >=B  < I >  THEN  7050 
IF  < A >X )  OR  <X>  =  X0>  THEN  7035 
T  a=A  + A 
M= I  NT ( Ta) 

IF  TaOM  THEN  7035 
J  = I  NT  < . 5*M> 

IF  A=J  THEN  7510 
GOTO  7535 
T 1 =A*LOG  <  X ) -X 
R=EXP<Tl)*FNGammai <A> 

GOTO  7115 
L=X/A 

IF  L=0  THEN  8025 
S= . 5-L+ . 5 
Z=FNR1 og(L) 

IF  2>=228/A  THEN  8015 
Y  =  A*Z 
R 1 =SQR  <  A ) 

IF  ABS ( S ) <  =E0/R 1  THEN  7830 
IF  ABS < S ) <  = . 4  THEN  7700 
T= 1 / <  A*A ) 

T 1  =  <  <  <  . 75*T-1 )*T  +  3. 5)*T-105)/<A*1260) 

T 1 =T 1 -Y 

R  =  Rpi *R1*EXP(T1) 

IF  R=0  THEN  8020 
IF  X<  =MAX <  A ,  Lrt  1 0 )  THEN  7140 
IF  X<X0  THEN  7605 
GOTO  7255 

******taylor  SERIES  FOR  P/R***** 

Apn=A+ 1 
T =X/Apn 
W< 1 >=T 

FOR  N=2  TO  20 
Apn=Apn+ 1 
T  =  T*  <  X^Apn ) 

IF  T<  *  1 E-3  THEN  7190 
W  <  N ) =T 
NEXT  N 
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N=23 
S  =  T 

To  1 = . 5*Ac  c 
Apn=Apn+ 1 
T=T*  <  X^Apn  > 

S  =  S  +  T 

IF  T >To 1  THEN  7200 
FOR  M 1 =N- 1  TO  1  STEP  -1 
S=S+W<M1 > 

NEXT  Ml 

flns  =  R/fl*  (.  1+S) 

Qans=. 5-fins+ . 5 
SUBEXIT 

!  *******ASYMPTOTIC  EXPANSION****** 

A 1 =  A- 1 
T  =  A1/-X 
W< 1 >=T 

FOR  N=2  TO  20 
A 1 =A 1 - 1 
T  =  T *  <  A  1 /X ) 

IF  ABS(TX  =  lE-3  THEN  7305 
W  <  N  >  =  T 
NEXT  N 
N  =  20 
S  =  T 

IF  ABS  ( T  X  =Ac  c  THEN  7335 

A  1 =A 1 - 1 

T=T*<A1/X) 

S  =  S  +  T 
GOTO  7310 

FOR  M=N- 1  TO  1  STEP  -1 
S=S+W<M> 

NEXT  M 

Qans=R/X*< 1+S) 

Ans= . 5-Qans+ .  5 
SUBEXIT 

!  *******TAYLOR  SERIES  FOR  P(A,X)/X**A****** 
A 1  =3 
C  =  X 

S  =  XXA  +  3> 

Tol =  3*AccXA+l ) 

A  1  =  A 1  +  1 
C=-C*<X/A1 ) 

T  =  CXA  +  A1  ) 

S  =  S  +  T 

IF  ABS ( T  )  >To 1  THEN  7390 
J  =  fl*X*<  <S^6-.  5XA  +  2)  >*X+1/<A+1  >  ) 

Z=A*LOG  <  X ) 

H=FNGam 1 (A) 

G=  1 +H 

IF  X< . 25  THEN  7450 
IF  A<  X/2 . 59  THEN  7475 
GOTO  7455 

IF  2  >- .  13394  THEN  7475 
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7455 

7460 

7465 

7470 

7475 

7480 

7485 

7490 

7495 

7500 

7505 

7510 

7515 

7520 

7525 

7530 

7535 

7540 

7545 

7550 

7555 

7560 

7565 

7570 

7575 

7580 

7585 

7590 

7595 

7600 

7605 

7610 

7615 

7620 

7625 

7630 

7635 

7640 

7645 

7650 

7655 

7660 

7665 

7670 

7675 

7680 

7685 

7690 

7695 

7700 

7705 

7710 

7715 

7720 

7725 

7730 


MW 


W=EXP<2> 

fins=W*G*< .5-J+.5) 

Qans*. 5 -flns+. 5 

SUBEXIT 

L«FNRexp<Z> 

W=. 5+L+. 5 

Qarrs=<U*J-L*Z)*G-H 
IF  Gians <  0  THEN  8040 
fins* . 5-Qans+ . 5 
SUBEXIT 

♦♦♦♦♦♦FINITE  SUMS  FOR  Q  WHEN  A>=1  FIND  2*A  IS  RN  INTEGER.  ******* 
S*EXP ( -X ) 

T  =  S 
k>l 
O0 

GOTO  7560 
T 1 =SQR  <  X ) 

S*FNErfc 1 <  T 1 . 0) 

T«EXP<-X)/<Spi *T1 > 

K  =  0 
C  =  - .  5 

FOR  N  =  K  TO  J-l 
C  =  C+1 
T=X*T/C 
S*S  +  T 
NEXT  N 
Qans=S 

fins*  .  5-Qans+  .  5 
SUBEXIT 

♦  ♦♦♦♦♦♦♦CONTINUED  FRACTION  EXPflNSION****^*^ 

Tol =MAX(5*E9, flee ) 

A2  1  =  1 
A2=  1 
B21-X 

B2=X+. 5-A+. 5 
C=1 

R2 1 =X*A2+C*A2 1 
B21=X^B2+C*B21 
Rmo*R21/B21 
C  =  C+1 
C 1 *C-R 

R2*R2 1 +C 1  * A2 
B2*B2 1 +C 1 *B2 
Ano*R2',B2 

IF  RBS(Rno-flfno)  >  =  To  1  *Rno  THEN  7635 

Qans=R*flno 

Ans*.5-Qans+.5 

SUBEXIT 

♦  ♦♦♦♦♦♦GENERAL  TEMME  EXPNS 

IF  <  ABS<SX=2*E9)  AND  C  A^E9^E9  >3 . 28E- 3  >  THEN  8055 
C=EXP(-V ) 

W=.5*FNErfcl(SGR<Y>,  1)  !♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦♦ 

U-l  /A 

2“SQR  <  2+  2  ) 

IF  L< 1  THEN  Z*-Z 
IF  I>*2  THEN  7785 


r» 

m 

I 


WW5 
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7735  IF  ABS <  S  >  <  =  1 E-3  THEN  7865 

7740  C0=<  CC<<<<<<<<  <D<  13>*Z+D<  12)  >*Z+D<  1 1  >  )  *Z+D  (1  0  )  )  *Z  +  D  < 9 >  )  *Z  +  D  (  8  >  >  *Z  +  D  <  7  >  >  * 

Z  +  D  <  6 ) >*Z+D<5> )*Z+D<4> )*Z+D<3) ) *Z  +  D < 2 ) ) *Z  +  D < 1 ) >*Z-Th 

7745  Cl  =  <  <<<<<<<<<  <D1 < 12>*Z  +  D1 < 1 1 >  >*Z  +  D1 < 10) >*Z  +  D1 < 9 > ) *Z  +  D 1 < 8 ) ) *Z  +  D 1 < 7 ) ) *Z  +  D 1 

<6> >  *Z  +  D 1 <5) ) *Z+D 1 <  4  ) ) *Z+D 1 <  3 ) )*Z+D1<2))*Z+D1(1))*Z+D10 

7750  C2=  <  < <  <  C  <  <  <  <D2< 10)*Z+D2<9) )*Z+D2<8) )*Z+D2<7> )*2+D2<6> )*Z+B2<5>  )*Z+D2<4) ) 

*2+D2<3) )*2+D2<2> ) *Z+D2  < 1 ) )*Z+D20 

7755  C3=( <  <  <  <  <  <D3<8)*Z+D3<7) )*Z+D3<6)  >*Z+D3(5> >*Z+D3(4) )*Z+D3<3) )*2+D3(2) >*2  + 

D3  <  1  )  > *Z+D30 

7760  C4= ( ( ( ( <D4<6)*Z+D4<5) )*Z  +  Il4<4> )*Z+D4<3> )*Z+D4<2)  )*Z+D4< 1 > )*Z+D40 

7765  C5=( <  <D5<4)*Z+D4<3) )*Z+D5<2) ) *Z+D5 ( 1 ) ) *Z+D50 

7770  C6=(D6<2)*Z+D6< 1 > >*Z+D60 

7775  T=<  <  <  <  <  <D70*U+C6)*U+C5)*U+C4>*U+C3)*U+C2>*U+C1 >*U+C0 

7780  GOTO  7940 

7785  IF  I >2  THEN  7815 

7790  C0=<  <  <  <  <D<6)*Z+B(5) )*Z+D<4) )*Z+D<3) >*Z+D<2) )*Z+D( 1 ) )*Z-Th 

7795  Cl  =  <  <  <D1 <4>*Z  +  B1 <3>  >*Z  +  D1 <2>  >*Z  +  D1 < 1 ) )*Z  +  D10 

7800  C2  =  D2  < 1 ) *Z  +  D20 

7805  T=<C2*U+C1 )*U+C0 

7810  GOTO  7940 

7815  T=<  <B<3)*Z  +  B<2) )*Z  +  B< 1 ) )*Z-Th 

7820  GOTO  7940 

7825  !  ***********TEMME  EXPANSION  FOR  L= 1  *********** 

7830  IF  A*E9*E9>3. 28E-3  THEN  8055 

7835  C=. 5-V+. 5 

7840  W=< . 5-SQR<Y)*<  . 5+< . 5-Y/3) )'Spi >/C 

7845  U= 1 / A 

7850  Z  =  SQR  <  Z  +  Z ) 

7855  IF  L< 1  THEN  Z=-Z 

7860  IF  I >=2  THEN  7910 

7865  C0=<  <  <  < ( <  B  ( 7 )*Z+B<6) >*Z+B(5> >*Z+B<4) )*Z+B(3> )*Z+B(2> )*Z  +  D( 1 ) )*Z-Th 

7870  C1=<<<<(B1(6)*Z+B1<5))*Z+B1<4>)*Z+B1(3))*Z+B1<2>)*Z+B1<1)>*Z+B10 

7875  C2=<  <  <  <B2(5)*Z+B2<4) )*Z+B2<3) >*Z+B2(2) >  *Z  +  B2 ( 1 ) >*Z+B20 

7880  C3=<  < (B3<4)*Z+B3<3) )*Z  +  D3<2) ) *Z+B3 ( 1 ) )*Z  +  D30 

7885  C4=<B4<2)*Z+B4< 1 >  >*Z+D40 

7890  C5=<B5<2)*Z+B5< 1 ) )*Z+B50 

7895  C6=D6< 1 )*Z+D60 

7900  T=( < ( <  <  <B70*U+C6)*U+C5)*U+C4)*U+C3)*U+C2)*U+C1 )*U  +  C0 

7905  GOTO  7940 

7910  IF  I >2  THEN  7935 

7915  C0=<B<2)*Z+B< 1 ) )*Z+C0 

7920  C1=B1 < 1 )*Z+B10 

7925  T=(D20*U+C1 )*U+C0 

7930  GOTO  7940 

7935  T  =  B  <  1  ) *Z-Th 

7940  IF  L< 1  THEN  7960 

7945  Qans  =  C*(W  +  Rpi *T/R1  ) 

7950  Ans= . 5-Qans+ . 5 

7955  SUBEXIT 

7960  Ans  =  C*<W-Rpi  ♦T.-'Rl  ) 

7965  Qans= . 5-Ans+ . 5 

7970  SUBEXIT 
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7975  !  ********«SPEC IRL  CflSES************ 

7980  IF  X >= . 25  THEN  8000 

7985  Rns=FNErf<SQR<X) ) 

7990  Qans=.  5-Flns  +  .  5 

7995  SUBEXIT 

8000  Qans=FNErfc 1 <SQR<X) , 0) 

8005  Rns=. 5-Qans+. 5 

8010  SUBEXIT 

8015  IF  RBS < S ) <  =2*E9  THEN  8055 

8020  IF  X>fl  THEN  8040 

8025  flns  =  0 

8030  Qans=l 

8035  SUBEXIT 

8040  Rns  = 1 

8045  Qans=0 

8050  SUBEXIT 

8055  flns=2 

8060  SUBEND 

8065  !  ********************RLNREL************************* 

8070  DEF  FNfll nrel <fl) 

8075  P 1 =- 1 . 29418923022E0 

8080  P2=. 405303492862 

8085  P3  =  -1.788?45*e>012E-2 

8090  Q1=-1.62752256355E0 

8095  Q2=. 74781 1014038 

8100  Q3=-8. 45104217946E-2 

8105  IF  RBS  <  R ) > . 375  THEN  8130 

8110  T=R/<fl+2) 

8115  T2=T*T 

8120  W=<  <  <P3*T2+P2)*T2+P1 )*T2+l)/< ( < Q3*T2+Q2 > *T2+Q 1 >*T2+1 > 

8125  RETURN  2*T*W 

8130  X= . 5+fl* . 5 

8135  RETURN  LOG<X) 

8140  FNEND 

8145  !  *******************GRMMR I  ************************************ 

8150  !  Gammai =1/<GRMMR  FUNCTION  OF  R).  !  NEEDS  GRM 1 . 

8155  DEF  FNGammai <fi) 

8160  IF  R >  =  20  THEN  8245 

8165  E  =  5E- 1 3 

8170  I  =  I  NT ( R ) 

8175  1 1 =fl- 1 

8180  Gammai=l 

8185  IF  I1O0  THEN  8200 

8190  16=0 

8195  GOTO  8220 

8200  I6=-.5+Il-.5 

8205  H=FNGam 1(11)  ! 16=11-. 5- - 5 

8210  Gammai = 1 1  +  1 1 *H 

8215  IF  1=0  THEN  RETURN  Gammai 

8220  12=1-1 

8225  16=16*1 

8230  Gammai =Gammai / 16 

8235  IF  I6CI2  THEN  8225 

8240  RETURN  Gammai 


*-i  «v*h  a*.  «.» v«  s*  vt»  scrran 
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8245 

8250 

8255 

8260 

8265 

8270 

8275 

8280 

8285 

8290 

8295 

8300 

8305 

8310 

8315 

8320 

8325 

8330 

8335 

8340 

8345 

8350 

8355 

8360 

8365 

8370 

8375 

8380 

8385 

8390 

8395 

8400 

8405 

8410 

8415 

8420 

8425 

8430 

8435 

8440 

8445 

8450 

8455 

8460 

8465 

8470 

8475 

8480 

8485 

8490 

8495 


IF  fl<  69 . 5  THEN  8260 

Gammai =0 

RETURN  Gammai 

B=. 918938533205 

Il=l/A 

12=11*11 

I6=<  .  5-A>*L0G(A>+A-B-< ( 1/1260*12-1/360) *12+1/1 2 >*  1 1 
Gammai =  EXP  < 1 6 ) 

RETURN  Gammai 
FNENIi 

****************RCOMP*************** 

COMPUTATION  OF  X'A*EXP < -X ) /GAMMA t A ) 


SUBROUTINES  NEEDED:  GAM  1 ,  GAMMAI,  RLOG. 

DEF  FNRcompCA, X) 

Rt2pi n=. 398942280401  ! 1/SQRC2*PI ) 

Rc  omp  =  0 

IF  A >  =  20  THEN  8360 
T  =  A*LOG ( X )  -X 
IF  A >= 1  THEN  8350 
Rcomp=A*EXP<T)*< 1+FNGaml <A> ) 

RETURN  Rcomp 

Rcomp=EXP<T)*FNGammai (A) 

RETURN  Rcomp 
U  =  X'fi 

IF  U=0  THEN  RETURN  Rcomp 
T* 1 / < A*A ) 

Tl  =  (  C  <  .  75*T-1 >  *T  +  3 . 5)*T-105>/< 1 260*A  > 

Tl=Tl-A*FNRlog<U) 

Rcomp=Rt2pin*SQR<A)*EXP(Tl> 

RETURN  Rcomp 
FNEND 

!  ****************C0MPUTES  Gamln<My)=LOG(GAMMA(My)**************** 
!  NEEDS  GAMLN1 
DEF  FNGam 1 n  <  My ) 

J5=. 418938533205  ! J5= . 5*L0G < 2*P I -  1  > 

J0=8. 33333333333E-2 
Jl=-2. 7777777777E-3 
J2=7. 93650663 184E-4 
J3=-5. 95156336429E-4 
J4=8. 20756370354E-4 
IF  My > . 8  THEN  8460  !Lg<My) 

Garni n=FNGaml nl <My)-L0G<My> ! SUBROUTINE  FOR  LOG < GAMMA < My+ 1 >  > 

RETURN  Gamin 
IF  My  >  =  2 . 25  THEN  8480 
Ny=- . 5+My- . 5 
Gamln=FNGamlnl<Ny) 

RETURN  Gamin 

IF  My >= 1 5  THEN  8535 

Ny= I  NT ( My- 1 . 25) 

Ay  =  My 
Wy  =  1 


m 
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8500  FOR  I  y=  1  TO  Hy 
8505  Ay = Ay- 1 

8510  Wy=Ay*Wy 

8515  NEXT  Iy 

8520  flyy=fly-l 

8525  Gam  1 n=FNGam 1 n 1 < flyy >  +  LOG C Wy > 

8530  RETURN  Gamin 

8535  Ayy=l/<My*My) 

8540  Wy=<  <  < ( J4*flyy+J3)*Ayy+J2>*Ayy+Jl >*flyy+J0)/My 

8545  Garni n=J5+Wy+(My-.5)*CL0GCMy)-l) 

8550  RETURN  Gamin 

8555  FNEND 

8560  !  *************COMF'UTES  R 1 og=L- 1 +LOG C L )  FOR  L>0. ****************** 

8565  DEF  FNRlog(L) 

8570  B3=. 333333333333 

8575  F2= .  2 

8580  F3=l . 42857142857E-1 

8585  F4«l . 1 1 1 1 1 1 11 11  IE-1 

b  5  :*  0  r5-3.G3o30  3  t)90^1L-2 

8595  F6=5. 88434334 1 73E-2  ! **F6  = . 695- 1 -LOG < . 695 > . ** 

8600  F7=4. 56512608816E-2  (**F7=4/3- 1 -LOG C 4/3 > . ** 

8605  S 1  = . 5-L+. 5 

8610  IF  <  L<  . 548 )  OR  <L>1.681)  THEN  8685  ! **SUBROUT I NE  FOR  R2  =  L- 1 -LOG C L > . ** 

8615  IF  L<. 83667  THEN  8640 

8620  IF  L  > 1 .  1547  THEN  8655 

8625  R2=-S 1 

8630  W1=0 

8635  GOTO  8665 

8640  R2=CL-. 695)/. 695 

8645  W1=F6-R2*. 305 

8650  GOTO  8665 

8655  R2=- . 75*S 1 - . 25 

8660  W 1 =F7+R2*B3 

8665  R2=R2/CR2+2> 

8670  T=R2*R2 

8675  Rlog=2*T*<l/C.5-R2+.5>-R2*<<<<F5*T+F4)*T+F3>*T*F2>*T+B3>>+Wl 

8680  RETURN  Rlog 

8685  R1 og=-Sl-LOGCL) 

8690  RETURN  Rlog 

8695  FNEND 

8700  !  **************GAM1=1/GAMMACA+1 >-l ,  -. 5< =A< = 1 . 5 . *************** 

8705  DEF  FNGam 1  <  A ) 

8710  DIM  P<7), Q<5), R(9) 

8715  PC  1  )  =  . 577215664902 

8720  P(2)=-. 409078193006 

8725  P(3)=-. 230975380858 

8730  P(4)=5. 97275330452E-2 

8735  P(5)=7. 66968181649E-3 

8740  P<6)=-5. 14889771324E-3 

8745  PC7)=5. 8959742861  IE-4 

8750  G  < 1 )  =  1 

8755  QC2>  =  . 427569613095 

8760  Q(3)=. 158451672430 

8765  G(4)=2.61 132021441E-2 
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8770  Q<5)=4. 23244297897E-3 

8775  R< 1 >=-. 422784335098 

8780  R<2>=-. 771330383816 

8785  R<3)=-. 244757765222 

8790  R<4)=. 1 18378989872 

8795  R<5>=9. 30357293360E-4 

8800  R(6)=-l . 18290993445E-2 

8805  R(7>=2. 23047661 158E-3 

8810  R<8)=2. 66505979059E-4 

8815  R<9>=-1 . 32674909766E-4 

8820  Sl=. 273076135304 

8825  S2=5. 59398236957E-2 

8830  T  =  fl 

8835  D  =  fi- . 5 

8840  IF  D  >0  THEN  T=B-.5 

8845  IF  T<  >0  THEN  8860 

8850  Gaml=0 

8855  RETURN  Garni 

8860  IF  T< 0  THEN  8895 

8865  Gam  1  =  < ( <  <  <  < P < 7 ) *T  +  P < 6 ) > *T  +  P < 5 > ) *T+P < 4 ) >*T  +  P<3>  >*T+P<2> )*T  +  P( 1 ))/(<( CQ(5> 

*T  +  Q<4) )*T  +  Q<3) )*T  +  Q<2) )*T  +  Q<  1  )  ) 

8870  IF  D >0  THEN  8885 

8875  Gaml=R*Gaml 

8880  RETURN  Garni 

8885  Gam  1 =T/fl* < - . 5  +  Gam 1 - . 5 ) 

8890  RETURN  Garni 

8895  Gam  1  =  <  < <  <  <<  <  < R < 9 ) *T  +  R < 8 > ) *T  +  R C 7 > ) *T  +  R ( 6 ) ) *T  +  R ( 5 ) > *T  +  R ( 4 > ) *1 +R ( 3 ) > *T  +  R ( 2 ) 

>*T  +  R< 1 ) )/<  <S2*T  +  S1)*T  +  1  ) 

8900  IF  L >0  THEN  8915 

8905  Gam  1 < . 5  +  Gam 1  + . 5 ) 

8910  RETURN  Garni 

8915  Gam  1  =T*Gam  1  /Fi 

8920  RETURN  Garni 

8925  FNEND 

8930  !  ♦**********4HHHf***REXP=<EXP<X)-l>/X***************** 

8935  DEF  FNRexp<X> 

8940  Pl=9. 1404191482E-10 

8945  P2=2. 3808236104E-2 

8950  Ql=-. 499999999086 

8955  Q2=. 107141568981 

8960  Q3=-l . 19041 179761E-2 

8965  Q4=5. 9513081 1 186E-4 

8970  IF  RBS ( X ) > . 15  THEN  8985 

8975  Rexp=< (P2*X  +  P1 )*X+1 )/<  <  <  < Q4*X  +  Q3 ) «X  +  Q2 ) *X  +  Q 1 > *X+ 1  ) 

8980  RETURN  Rexp 

8985  W  =  EXP  <  X ) 

8990  IF  X>0  THEN  9005 

8995  Rexp=<W-. 5-. 5)/X 

9000  RETURN  Rexp 

9005  Rexp=W*<.5+(.5-l/W))/X 

9010  RETURN  Rexp 

9015  FNEND 
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9020  !  **********************ERF<X)********************* 

9025  DEF  FNErf(X) 

9030  X2=X*X 

9035  IF  flBS <  X  )  > .  5  THEN  9050 

9040  Erf=X*( < (-3. 560984370 18E-2*X2+6. 99638348862 >*X2  +  21 . 97926 1 6 1 83 ) *X2+242 . 66 
7955231 )/<< (X2+ 15. 0827976304 )*X2+91 . 1 649054045 > *X2+2 1 5 . 05887587 > 

9045  RETURN  Erf 

9050  K4  =  flBS  <  X ) 

9055  IF  K4 >4  THEN  9085 

9060  K3=( <  < (-6. 0858 151 9597E-6*K4+. 56 4 371 606 864) *K4+4. 267720 10709) *K4+1 4. 57189 

85969 )*K4+26. 0947469561 )*K4 +22. 8989928517 

9065  K 3  =  K 3  -'  <  CCU K4+7 . 56884822936 ) *K4+26 . 2887957588 ) *K4+50 . 2732028638 >  *K4  +  5 1 . 9 

335706876 )*K4+22. 8989857499) 

9070  Erf=l-EXP(-X2)*K3 

9075  IF  X<0  THEN  RETURN  -Erf 

9080  RETURN  Erf 

9085  Er f  = 1 

9090  IF  K4>=4. 883  THEN  9115 

9095  K6=1'X2 

9100  Bl  =  . 564189583548  !  1/SQRCPI) 

9105  Erf=<Bl-K6*< < ( 3 . 243 1 95 1 9278E-2*K6+ . 2439 1 1 029489 ) *K6  + . 1 1 9903955268 ) *K6+ . 0 

12130827639 )✓<  < (K6+1 . 4377 1 227937 ) *K6+ . 48955244 1 96 1 ) *K6  +  4 . 30026643453E-2 ) ) /K4 

9110  Erf=l-EXP<-X2)*Erf 

9115  IF  X<0  THEN  RETURN  -Erf 

9120  RETURN  Erf 

9125  FNEND 

9130  !  **********»**»*»*ERFC1 <X>******************* 

9135  !  IF  I nd=0  THEN  Erf c 1 =ERFC < X ) ;  ELSE  Er f c 1 =EXP < X*X > *ERFC < X ) . 

9140  DEF  FNErfc 1 (X, Ind) 

9145  X2=X*X 

9150  IF  flBS<X)>.5  THEN  9170 

9155  Erf  c 1  =  . 5-X*  C (( -3 . 560984  370 1 8E-2*X2  +  6 . 99638348862 ) *X2  +  2 1 . 9792616183>*X2  +  2 

42.667955231  )/( ( ( X2+ 1 5 . 0827976304 ) *X2  +  9 1 .  1 64  905404  5 ) *X2  +  2 1 5 . 05887587 )  +  .  5 

9160  IF  Ir>dO0  THEN  RETURN  EXP  <  X2  )  *Er  f  c  1 

9165  RETURN  Erf c  1 

9170  K4=RBS<X) 

9175  IF  K4  >4  THEN  9205 

9180  K3=<  < ( (-6. 08581519597E-6*K4+. 56437 1606864 >*K4  +  4. 267720 1 0709 )*K4+14. 57189 

85969 )*K 4+26. 0947469561 ) *K4+22 . 89899285 1 7 

9185  Erf c l=K3/< (<( (K4+ 7. 56884822936 ) *K4 +26 . 2887957588 > *K 4+50 . 2732028638 ) *K 4 +5 

1 . 9335706876>*K4+22. 8989857499) 

9190  IF  I nd=0  THEN  9245 

9195  IF  X<0  THEN  RETURN  2*EXP < X2  ) -Er f c  1 

9200  RETURm  Erfcl 

9205  IF  XC-4.89  THEN  9260 

9210  IF  (X>  =14. 9891301)  FIND  <lr>d  =  0)  THEN  RETURN  0 
9215  K6=l/X2 

9220  Bl  =  . 564189583548  1  1/SQR<PI) 

9225  Erf  c  1  =  <B1-K6*  (.(.<  3. 2431 951  9278E-2*K6+.  24391  1029489)  *K6+.  1  1990  3955268  )*K6  + 

. 012130827639) ( ( K6+1 . 4377 1 227937 ) *K6+ . 48955244 1 96 1 ) *K6+4 . 30026643453E-2 ) > 'K4 


9230 

IF  Ind 

=0  THEN  9245 

9235 

IF  X<-0 

THEN  RETURN  2*EXP ( X2 ) -Erf c  1 

9240 

RETURN 

Erfc  1 

9245 

Erfc  1  = 

EXP< -X2)*Erf c 1 

9250 

IF  X<0 

THEN  RETURN  2-Erfc 1 

9255 

RETURN 

Erfc  1 

C-  !0 


>0  THEN  RETUR 
2 

♦♦COMPUTES  Ga 

1 n 1 ( flyy  ) 

215664902  ! 

203922187 

8860593647 

0427615534 

2055799310 
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k4  =  3 7. 390276351  2.-  <  Daws +  4. 07068101667+K4) 

Daws=  <  .  214221965778.-'<Daws-4. 9 1 60536574 1 +K4 )+ . 500652754437 )/« 1 
RETURN  Daws 

IF  ABS  <  K  1  )  >  =  5  THEN  9545 

K4=  17. 76  1778077 ''<  Daw s-4.  108  3233770) 

K4  =  -.  537  1427298  l''<  Daws -10. 38502482  1+K4) 

Daws  =  <  .  2492463242 l/< Daws- 1 . 5853935006  +  K4 )  +  . 50000965 1 99 ) 'K 1 
RETURN  Daws 

K4=-4.  1  2544065608-'  <  Daws-  1  1  .  1952164237) 

K4  =  -2. 48787658804  - '(Daws-4. €731 20221 4  1 +K4) 

Daws*.  5*<l  +  <.  749999  1  90567.-'<  Daws-2. 508171  1 6686+K4 > + .  500000001674  >/Daw 

RETURN  Daws 
FNEND 

!  THIS  PROGRAM  IS  CALLED  "ELLCV" .  IT  SUPPLIES 
i  THE  ELLIPTICAL  COVERAGE  FUNCTION:  El  1 c  w  =  P ( R ,  H, K , SI , S2> . 

•  El  lew  DENOTES  THE  PROBABILITY  OF  A  SHOT,  NORMALLY  DISTRIBUTED  WITH 
!  MEAN  (0,0)  AND  STANDARD  DEVIATIONS  S1,S2  I1'  THE  X  AND  Y  DIRECTIONS, 

!  RESPECTIVELY,  FALLING  IN  A  CIRCLE  IN  THE  XY-PLANE  OF  RADIUS  R  AND 
i  CENTERED  AT  ( H , K ) .  THE  INPUT  IS  R,H,K,S1,S2.  THE  OUTPUT  IS  Ellcv. 
i  PROGRAM  IS  SET  FOR  8-DECIMAL-DIGIT  ACCURACY  IN  El  lew. 

!  ELLCV  USES  ERFC,  ERF  WITH  12  DIGIT  RELATIVE  ACCURACY. 

!  SOURCES:  NWL  REPORT  #1710,  AUG. I960.  MATH  OF  COMP.  OCT.  1961, 
i  PP.  375,382.  NSWC  REPORT  #83-13,  NOVEMBER,  1982. 

DEF  FNE 1 1 cw<R, H, K, SI , S2> 

COM  X<*),Y(*> 

Bl=. 564189583548 !  l/SQRCPI) 

Z3=. 000000O5*S1*S2 
IF  R*R<  =Z3  THEN  RETURN  0 
A=6 . 027  '6  '5.6123 

8=1.41421356237  !  SQRC2) 

B2  =  36  '31.49791129  ! B2=A*A 

h2=H*H+K*K 

H8=ABS ( H ) 

K  8  =  A  B  S  (  K  ) 

D=MAX( SI , S2 ) 

T  =  R- A*D 

'  ***PP0CEED  TO  SEE  IF  P=0  OR  P=1 
IF  T  < 0  THEN  9695 
IF  T  *  T  <  H  2  THEN  9695 
RETURN  1 

IF  R-H8  +  A*S 1 <  =0  THEN  RETURN  0 
IF  R-K8+A*S2<  =0  THEN  RETURN  0 
S0=SQR ( H2 ) 

IF  S0<  =R  THEN  9730 
D=  <S0-R)*(S0-R>^(D*D ) 

IF  R*R*EXPC- . 5+D ■ >Z3  THEN  9730 
RETURN  0 

IF  S10S2  THEN  9750 
H8  =  S0 
K8  =  0 

IF  R-H8+A*S1  =0  THEN  RETURN  0 
IF  H8*K3=0  THEN  9855 


RETURN 

RETURN 


RETURN 


M  '  'll  I 


9755 
9760 
9765 
9770 
9775 
9780 
9785 
9790 
9795 
9800 
9805 
9810 
9815 
9820 
9825 
9830 
9835 
9840 
9845 
9850 
9855 
9860 
9865 
9870 
9875 
9880 
9885 
9890 
9895 
9900 
9905 
9910 
9915 
9920 
9925 
9930 
9935 
9940 
9945 
9950 
9955 
9960 
9965 
9970 
9975 
9980 
9985 
9990 
9995 
10000 
10005 
10010 
10015 
10020 
)  +  5  *  S  9  ) 


RETURN  0 
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H9=H8/S1 

K9=K8/S2 

D=H9*H9+K9*K9 

Z2=R/S2 

IF  D<  =B2  THEN  9855 
Z 1 =R+ A*M I N ( S 1 , S2> 

IF  H2<  Z 1 *Z 1  THEN  9855 

Q=S2/S1 

Q1=Q*Q 

F=Q1»H9*H9+K9*K9 

Z 1 =Z2*Z2*F/D 

Z=D-Z 1 -B2 

IF  Z<0  THEN  9825 

IF  Z*Z-4*Z1*B2>=0  THEN  RETURN  0 

T1=H8*H8+Q1*K8*K8 

Z8=B2*S1*S1*T1/H2 

R2=R*R 

V=H2-R2-Z8 

IF  V<0  THEN  9855 

IF  Y*Y-4*R2*Z8>=0  THEN  RETURN  0 

28  =  0  !  ****F  I  Nil  LIMITS  OF  INTEGRATION 

2=K8+A*S2 

H3=K8-A*S2 

S0  =  S1 

S9  =  S2 

Z  =  R-Z 

H5  =  0 

D  1  =0 

IF  Z>  =  0  THEN  Dl=SGlR<Z/R) 

IF  H3>=0  THEN  9915 
E3= 1 -D 1 
GOTO  9925 
E3=SQR< 1-H3/R>-D1 
H5-1 

IF  Z8<  >0  THEN  9965 
Z8=  1 
F  =  E3 
T  =  D  1 

Z=H8+A*S1 
H6  =  H5 

H3=H8-A*S 1 
GOTO  9880 

IF  F  >=E3  THEN  10010 

E3  =  F 

D  1  =T 

S9  =  3  1 

Z8  =  H8 

S0  =  S2 

H8  =  K8 

K8  =  Z8 

H5  =  H6 

E3= . 5*E?  '♦♦♦BEGIN  GAUSSIAN  I NTEGR AT  I  ON* ♦ ♦ ♦ ♦ 

N  =  3 

IF  <  H8-R  >-3 . 5*S0  -i  OR  <  K8-R  >-3 . 5^S9  >  THEN  N=E3*R^  (  .  34/S0+  1  ✓  <  .  01 


INTEGRAT ION***** 


:'*ABS^R-K8 


C-i'i 


to 


10025 

10030 

10035 

10040 

10045 

10050 

10055 

10060 

10065 

10070 

10075 

10080 

10085 

10090 

10095 

10100 

10105 

10110 

10115 

10120 

10125 

10130 

10135 

10140 

10145 

10150 

10155 

10160 

10165 

10170 

10175 

10180 

10185 

10190 

10195 

10200 

10205 

10210 

10215 

10220 

10225 

10230 

10235 

10240 

10245 

10250 

10255 

10260 

10265 

10270 

10275 

10280 

10285 

10290 

10295 
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Z2=R/<B*S9> 

R8=R/ ( B*S0 ) 

H8=H8/CB*S6> 

K8=K8^<B*S9) 

IF  N<  2 . 75  THEN  10065 
J  =  3 1 
N  1  =  1  2 

GOTO  10155 
IF  N< 1 . 35  THEN  10085 
J  =  2  1 
N  1  =  1  0 

GOTO  10155 
IF  N< . 75  THEN  10105 
J=  1  3 
N  1  =8 

GOTO  10155 
IF  N< . 35  THEN  10125 
J  =  7 
N  1  =6 

GOTO  10155 
IF  N< . 15  THEN  10145 
J  =  3 
N  1  =4 

GOTO  10155 
J  =  0 
N  1  =3 
Z  =  0 

Y=B 1 *E3*R8 

K9=0  !  1 . 04E-8 

H9  =  2  11.9999999702 

G3  =  0 

IF  K8=0  THEN  10390 
23  =  0 

FOR  I =-N 1  TO  N 1 

IF  1=0  THEN  18370 

T=E3* ( SGN ( I >*XC J+RBS< I ) >+l >+Dl 

T9=T*T 

T1=R3*< 1-T9) 

T2=<H8-T1 >*<H8-T1 ) 

T4=EXP ( -T2 ) 

IF  H8<  >0  THEN  10240 

T4  =  T  4  +  T  4 

GOTO  10255 

IF  H5<  >0  THEN  10255 

T2=  <  H8  +  T 1 )*(H8+T1 ) 

T4=T4+EXP(-T2) 

IF  Z=0  THEN  10270 
K5  =  H9 

GOTO  10365 
Z1=Z2*T*SQR<2-T9) 

K1=K8-Z1 

IF  K  1  <  =  -R  THEN  10300 
Erfc=FNErfc  <  K 1 > 

K5=Erf c 
GOTO  10315 


II7P» 

m 


wJ!J 
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10300  2=1 

10305  K5  =  H9 

10310  GOTO  10365 

10315  IF  23=0  THEN  10330 

10320  K5=K5-K9 

10325  GOTO  10365 

10330  K 1 =K8+2 1 

10335  IF  K  1  >  =  fl  THEN  10355 

10340  Erfc=FNErfc <K1) 

10345  K5=K5-Erfc 

10350  GOTO  10365 

10355  K5=K5-K9 

10360  23=1 

10365  G3=G3+K5*T4*T*Y<J+RBS< I > > 

10370  NEXT  I 

10375  E 1 1 c  v=Y*G3 

10380  IF  El  1 c  v  > 1  THEN  RETURN  1 

10385  RETURN  El  1cm 

10390  FOR  I  =  —  N 1  TO  N1 

10395  IF  1=0  THEN  10370 

10400  T  =  E3*  ( SGN  <  I  )*X<  J+RBSCI  >  >  +  l  )  +  Dl 

10405  T9=T*T 

10410  T 1 =R8* ( 1-T9) 

10415  T2=<H8-T1 )*(H8-T1 ) 

10420  T 4  =  EXP  < -T2 ) 

10425  IF  H5<  >0  THEN  10440 

10430  T2=<H8  +  T1 )*<H8+T1  ) 

10435  T4  =  T4  +  EXPOT2) 

10440  IF  2=0  THEN  10455 

10445  K5=H9 

10450  GOTO  10365 

10455  K1=-22*T*SQR<2-T9) 

10460  IF  K 1  <  =-F)  THEN  10480 

10465  Erfc=FNErf <-Kl )  ! FNErf  c ( K 1 ) 

10470  K5=2*Erfc  !2*<Erfc-l> 

10475  GOTO  10365 

10480  2=1 

10485  K5=H9 

10490  GOTO  10365 

10495  FNEND 

10505  DEF  FNEr f  c  <  X  > 

10510  K6=X*X 

10515  IF  BBS  <  X ) > . 5  THEN  10525 

10520  RETURN  . 5-X* <  <  < -3 . 560984370 1 SE-2*K6  +  6 . 99638348862 > *K6  +  2 1 . 97926 1 6 1 83 > *K6  + 
242.667955231 >/<<< K6+ 1 5 . 0827976304 ) *K6+9 1 . 1 649054045 ) *K6  +  2 1 5 . 05887587 >  +  . 5 
10525  K4=RBS ( X ) 

10530  IF  K4  >4  THEN  10550 

10535  Erf c  =  <  <<  <-6. 08581519597E-6*K4  +  . 564371 606864 >*K4  +  4. 267720 1 0709 > *K4  + 1 4 . 571 

8985969 )*K4+26. 0947469561 ) *K4+22 . 89899285 1 7 

10540  Erfc=EXP<-K6)*Erfc  ••'<<<<  (K4  +  7. 56884822936)  *K4 +  26. 2887957588  >*K4  +  50. 273202 

8638 ) *K4+5 1 . 9335706876 > *K4 +22 . 8989857499 > 


vS-*i 

y. 
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10545  GOTO  10580 

10550  IF  XC-4.89  THEN  RETURN  2 

10555  IF  X>=14. 9891301  THEN  RETURN  0 
10560  K6=1^K 6 

10565  Bl=. 564189583548 !  1/SQR(PI> 

10570  Erfc=<Bl-K6*( <<3.24319519278E-2*K6+.24391 1 029489 ) *K6+ . 1 1 9903955268 ) *K6+ . 

01 21 30827639 )/< ( CK6+1 . 4  377 1 227937 ) *K6  + . 48955244 1 96 1  > *K6  +  4 . 300266434  53E-2 >  > /K4 

10575  Erfc=EXP(-X*X)*Erfc 

10580  IF  X<0  THEN  RETURN  2-Erfc 

10585  RETURN  Erfc 

10590  FNEND 

10595  !  *****************CIRCV****************** 

10600  !  THIS  SUBROUTINE  IS  CALLED  CIRCV.  INPUT:  R,D,S,V. 

10605  !  IF  V#0.  OUTPUT:  C i rc v=C I RCULAR  COVERAGE  FUNCTION.  Circv  GIVES  THE 

10606  !  PROBABILITY  OF  A  SHOT  FALLING,  UNDER  A  NORMAL  DISTRIBUTION  WITH 

10610  !  MEAN  (0,0)  AND  EGUAL  STANDARD  DEVIATIONS,  S,  IN  A  CIRCLE  OF  RADIUS  R 

10611  !  OFFSET  A  DISTANCE  D  FROM  (0,0). 

10615  !  IF  V=0 ,  OUTPUT:  C l rc v  =  GENERAL I  ZED  CIRCULAR  ERROR  FUNCTON.  Circu  GIVES 

10616  !  THE  PROBABLITY  OF  A  SHOT  FALLING,  UNDER  A  NORMAL  DISTRIBUTION  WITH 

10620  !  MEAN  (0,0)  AND  STANDARD  DEVIATIONS  SI  AND  S,  IN  A  CIRCLE  OF  RADIUS  R 

10621  !  CENTERED  AT  (0,0).  D  =  Sl/S  <=  1. 

10625  '  FOR  V  =  0,  WITH  SI  =  0,  S1#S,  C i rc u=ERF < R/ < SQR < 2 > *S ))  . 

10630  !  SUBROUTINES  NEEDED:  ERFC2 ,  ERF. 

10635  !  THE  PAR  I AL  DERIVATIVE  OF  Ci rev  WITH  RESPECT  TO  R  IS  CONTAINED  WITHIN  THE 

10636  !  ROUTINE  AND  CAN  BE  EXTRACTED  IF  DESIRED. 

10640  !  SOURCES:  MATH  OF  COMP  APRIL  1 96 1 , PP 1 69 , 1 73  AND  OCT. 1961,  PP  375,382. 

10641  !  NWL  REPORT  #1768,  JAN.  1962.  NSWC  REP0RT#83- 1 3 ,  NOV.  1982. 

10645  !  IEEE  TRANS.  INFO.  TH.  APRIL  1965,  P.  312. 

10650  !  PROGRAM  IS  SET  FOR  APPROXIMATELY  EIGHT  DECIMAL  DIGIT  ACCURACY. 

10655  DEF  FNCi rcvCR, D, S, V) 

10660  IF  R=0  THEN  RETURN  0 
10665  Bl=. 564189583548 !  l/SQRCPI) 

10670  C7=. 707106781 187 

10675  C 1 = 1 E-8 

10680  R 1  =R/S 

10685  IF  V< >0  THEN  10740 

10690  IF  ABS < D- . 5 ) <  = . 5  THEN  10700 

10695  PRINT  "ERROR:  D>1  OR  D<0,  IN  CIRCV." 

1 C~00  D1=P 

10705  R1=R/S 

10710  IF  D 1 <  >0  THEN  10720 

10715  RETURN  FNEr f  <R1*C7>  !  V  =  0  AND  D  =  S1=0. 

10720  T=.5*R1/D1 

10725  R 1 =T*  < 1+D1 ) 

10730  D1=T*< 1  —  D 1 ) 

10735  GOTO  10745 
10740  D1=D/S 
10745  A  1 =R 1 -D 1 
10750  K 1 = ABS  <  A 1 ) 

10755  IF  Kl<6. 02738  THEN  10770  15.6123 

10760  IF  A 1 <  0  THEN  RETURN  0 
10765  RETURN  1 


b 

i 
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v 


NSWC  TR  87-27 


10770 

10775 

10780 

10785 

10790 

10795 

10800 

10805 

10810 

10815 

10820 

10825 

10830 

10835 

10840 

10845 

10850 

10855 

10860 

10865 

10870 

10875 

10880 

10885 

10890 

10895 

10900 

10905 

10910 

10915 

10920 

10925 

10930 

10935 

10940 

10945 

10950 

10955 

10960 

10965 

10970 

10975 

10980 

10985 

10990 

10995 

11000 

11005 

11010 

11015 

11020 

11025 

11030 

11035 

11040 


T=R 1 *D 1 
T3= . 5*R 1 *R 1 
B=. 5*D1*D1 
N  =  0 

IF  T >7  THEN  10930 
T 1 =C7*T- 1 
T2=T3*B 

50  =  EXP  < -T3-B ) 

S 1 =EXP  < -B  > 

IF  T3 > . 0005  THEN  10830 

S1=S1*T3*< . 5+<  <-T3/24+l/6)*T3-.  5>*T3+.  5  > 

GOTO  10835 

S1=S1-S0 

S2  =  S0 

T0  =  S  1 

N  =  N+ 1 

M=l/N 

S0=T2*M*M*S0 

T0=B*M*T0-S0 

5 1  =S 1 +T0 
S2=S2+S0 

IF  TON  THEN  10845 
IF  T0  >C 1  THEN  10895 
P  =  S1 

GOTO  11080 
N=N+  1 
M=l/N 

S0=T2*M*M*S0 
T0=B*M*T0-S0 
S 1 =S 1 +T0 
S2=S2+S0 
GOTO  10880 
T 1 =2*HBS  <  T3-B ) 

K 1 =K 1 *C7 
T3= 1 / <  T  +  T ) 

T2  =  SG!R<T3) 

Sl=. 5*fl 1 *R 1 

52  =  EXP  < -S 1  ) 

S0=B 1 *T2*S2 

T0=<R1+D1 >*C7*T2*FNErfc2<Kl , S2> 

T2=S 1 *T3 

T3=.5*T3 

S  1  =T0 

S2  =  S0 

N  =  N  +  2 

M  =  N- 1 

K1=M/N 

S0=K 1 *T3*S0 

T0=T1*S0-T2*K1*T0 

S0=M*S0 

S 1 =S 1 +T0 

S2=S2+S0 

IF  T0-C1>0  THEN  10990 
IF  S0-C1<=0  THEN  11075 
Tl«flBS<S0> 


C-37 


l. 1 **>  k*  i*'  J 
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11045 

11050 

11055 

11060 

11065 

11070 

11075 

11080 

11085 

11090 

11095 


N  =  N  +  2 
M  =  N-  1 

S0=M*M*T3*S0/N 

IF  ABS  <  S0  ) >  =  T 1  THEN  11075 

S2=S2+S0 

GOTO  11035 

P=. 5*ABS< 1+SGN(A1 )-S2-SGN(Al )*S1 ) 
IF  VO0  THEN  11090 
P=ABS(P+P+S2-1 ) 

IF  P>1  THEN  RETURN  1 
RETURN  P  ! ***EXI T*** 


11100  FNEND 

11105  !  ***************ERFC2**************** 

11110  !  Erfc=ERFC(X) .  Exp=EXPC-X»X)  IF  AVAILABLE,  OTHERWISE  Exp=-1. 

11115  DEF  FNErfc2(X, Exp) 

11120  K6=X*X 

11125  IF  ABS < X ) > . 5  THEN  11135 

11130  RETURN  . 5-X*< ( (-3. 5609843701 8E-2*K6+6. 99638348862 )*K6+21 . 97926 1 6 1 83 > *K6+ 
242. 667955231  <  <K6+ 15. 0827976304 ) *K6+9 1 .  1 64  9054045 ) *K6  +  2 1 5 . 05887587 )  +  . 5 

11135  K4  =  ABS  <  X ) 

11140  IF  K4  >4  THEN  11165 

11145  IF  Exp< 0  THEN  Exp=EXP(-K6) 

11150  Er f  c  =  <  < (<-6. 0858 151 9597E-6*K4+ . 56437 1 606864 ) *K4+4 . 267720 10709) *K4+ 14. 571 

8985969 )*K4+26. 0947469561 )*K4+22. 8989928517 

11155  Erfc=Exp*Erf'c/<  (  <  <  <  K4  +  7 . 56884822936  )  *K4  +  26 . 2887957588  )  *K4  +  50 . 27  32028638  ) 

*K4+51 .9335706876 )*K4+22. 8989857499) 

11160  GOTO  11200 

11165  IF  X< -4 . 89  THEN  RETURN  2 

11170  IF  X>=14. 9891301  THEN  RETURN  0 
11175  IF  Exp=- 1  THEN  Exp=EXP(-K6) 

11180  K6=l/K6 

11185  Bl=. 564189583548 !  lxSQR<PI) 

11190  Erfc=(Bl-K6*<  <  <3. 2431 951 9278E-2*K6+. 24391 1 029489 )*K6+. 1  1 9903955268 ) *K6+ . 

01 21 30827639) /<  <  <K6*1 . 4377 1 227937 ) *K6+ . 48955244 1 96 1  )*K6  +  4. 30026643453E-2 ) ) xK4 

11195  Erf  c  =Exp*Er f  c 

11200  IF  X<0  THEN  RETURN  2-Erf c 

11205  RETURN  Erfc 

11210  FNENB 


11215 

11220 

11221 

11225 

11230 

11235 

11236 
11240 

11245 

11246 

11250 

11251 

11255 

11256 
11260 


***************gquad****************** 

Gquad  IS  A  QUADRATURE  ROUTINE  WHICH  USES  GAUSSIAN  MULTIPLIERS  TO  OBTAIN 
THE  INTEGRAL  I  FROM  A  TO  B  OF  Fn(T),  (SEE  NEXT  SUBPROGRAM.) 

A , B  :  ARE  LOWER  AND  UPPER  LIMITS  OF  INTEGRATION. 

I  err  :  IF  0  THEN  ANSWER  O.K.  IF  1  THEN  MAXIMUM  NO.  OF  SUBDIVISIONS  USED. 
M  :  MAXIMUM  NO.  OF  EQUAL  SUBDIVISIONS  OF  [  A  ,  B  3  DESIRED  WITH  GAUSSIAN 
APPLIED  ON  EACH  SUBDIVISION. 

II  :  VALUE  OF  INTEGRAL  STORED  IN  II. 

E  :  |  I  1  —  1 0  |  <  E  SAYS  DIFFERENCE  BETWEEN  TWO  CONSECUTIVE  CALULAT I ONS  FOR 

I  ARE  LESS  THAN  E. 

N  :  2N  =  0RDER  OF  GAUSSIAN  MULTIPLIERS  USED.  ONE  MAY  SPECIFY  N=3,4,6,8, 
10,12.  ONLY  THESE  VALUES  OF  N  ALLOWED. 

X<*),Y<*)  :  STORED  VALUES  OF  GAUSSIAN  ABSCISSAS  AND  MULTIPLIERS  ON 

C-1,13. 

J  :  NO.  OF  EQUAL  SUBDIVISIONS  USED.  J  <=M. 
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11265  !  NEEDS  X<*>,  Y<*>  ARRAYS  WHICH  ARE  STORED  IN  "GAUSS";  FUNCTION  FN. 
11270  SUB  GquadCA, B, Ierr, M, 1 1 , E , N , N2 , R3 , H , K , L , S 1 , S2 ,  S3 ,  V ,  D ,  W9 > 

11275  COM  X<*),Y(*) 

11280  Jl=lerr=0 

11285  IF  N=3  THEN  11355 

11290  J 1  =3 

11295  IF  N=4  THEN  11355 
11300  J 1 =7 

11305  IF  N=6  THEN  11355 
11310  J 1  =  1 3 

11315  IF  N=8  THEN  11355 
11320  J 1 =2 1 

11325  IF  N= 1 0  THEN  11355 
11330  J  1  =3 1 

11335  IF  N= 1 2  THEN  11355 

11340  PRINT  "N  AND  J1  NOT  RIGHT" 

11345  BEEP 

11350  STOP 

11355  D3  =  B- A 

11360  J  =  0 

11365  1 0= 1 E50 

11370  J=J+1 

11375  E  1  =  A 

11380  G  =  0 

11385  D1=D3/J 

11390  D2=Dl/2 

11395  K 1 =0 

11400  K 1 =K 1 + 1 

11405  E0=E 1 

11410  E 1 =E0+D 1 

11415  E2=<E1  +E0>*.  5 

11420  FOR  I = -N  TO  N 

11425  IF  1=0  THEN  11445 

11430  T=D2*(SGN(I)*X<J1+ABS<I)>>+E2 

11435  F  =  Y ( J 1 +ABS  <  I  ) ) *FNFn < T , R3 , H , K , L , S 1 , S2,S3,V,D,W9> 

11440  G=G+F 

11445  NEXT  I 

11450  IF  KIOJ  THEN  11400 

11455  I  1 =D2*G 

11460  IF  ABS <  1 1-I0XE  THEN  11480 

11465  10=11 

11470  IF  JOM  THEN  11370 
11475  I er r = 1 

11480  N2=N*J*< J+l ) ! NO.  OF  FCN  EVALUATIONS. 

11485  SUBEND 

11495  !  NEEDED  FOR  GQUAD .  NEEDS  ELLCV,  CIRCV. 

11500  DEF  FNFn<T5, R3, H, K, L, SI , S2, S3, V, D, W9> 


11505 

COM  XC*> , Y<*> 

11510 

C8=l . 41421356237 

11515 

J  1  =0 

11520 

P=L-C8*S3*T5 

11525 

J 1 =EXP( -T5*T5 ) 

11530 

IF  L<  >0  THEN  11545 

11535 

J1=J1+J1 
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11540  GOTO  11565 

11545  IF  149  0  0  THEN  11565 

11550  V2  =  2*L*P''(S3*S3> 

11555  IF  V2  >27  THEN  11565 

11560  J1=J1*< . 5+EXP<-V2>+. 5) 

11565  R=SQR<R3*R3-P*P) 

11570  IF  V  =  2  THEN  11585 
11575  P=FNCi rcw(R, D, SI , V) 

11580  GOTO  11590 

11585  P=FNE1 1cw(R,H,K,Sl,S2) 

11590  J1=J1*P 

11595  RETURN  J1 

11600  FNEND 

11605  !  ***************ERF3**************** 

11610  !  Erf=ERF<X).  IF  AVAILABLE  Exp=EXP  < -X*X  >  ELSE  Exp=-1. 

11615  DEF  FNErf3<X, Exp) 

11620  X2=X*X 

11625  IF  ABS ( X ) > . 5  THEN  11640 

11630  Erf  =  X*<  << -3. 560984370 18E-2*X2  +  6. 99638348862 )*X2  +  21 . 97926 1 6 1 83 ) *X2+242 . 66 

7955231  )x<  <  <X2  + 15. 0827976304 ) *X2  +  9 1 . 1 649054045 ) *X2  +  2 1 5 . 05887587 ) 

11635  RETURN  Erf 

11640  K4  =  ABS  <  X ) 

11645  IF  K4  >4  THEN  11680 

11650  K3= ( < < (-6. 08581519597E-6*K4+. 56437 1 606864 ) *K4+4 . 267720 10709 >*K4+14. 57189 

85969 )*K4+26. 0947469561 )*K4+22. 8989928517 

11655  K3  =  K3/( < < <  < K4  +  7 . 56884822936 ) *K4  +  26 . 2887957588 ) *K4  +  50 . 2732028638 ) *K4  +  5 1 . 9 

335706876 )*K4 +22. 8989857499) 

11660  IF  Exp<  0  THEN  Exp  =  EXP(-X2) 

11665  Erf = 1 -Exp*K3 

11670  IF  X<  0  THEN  RETURN  -Erf 

11675  RETURN  Erf 

11680  Erf = 1 

11685  IF  K4  >  =  4 . 883  THEN  11715 
11690  K6=l/X2 

11695  Bl  =  . 564189583548  !  lxSQRCPl) 

11700  Erf=(Bl-K6*<  <  <3. 24319519278E-2  +  K6+. 24391  1 029489 ) *K6+ . 1 1 9903955268 ) *K6+ . 0 
1 21 30827639 )/<<  <K6+1 . 43771 227937)*K6+. 489552441961 ) +K6  +  4 . 30026643453E-2 ) ) /K4 
11705  IF  Exp<  0  THEN  Exp  =  EXP < -X2 ) 

11710  Erf*l-E\p*Erf 

11715  IF  X<0  THEN  RETURN  -Erf 

11720  RETURN  Erf 

11725  FNEND 

11730  !  ***************ELLRC**************** 

11735  !  FINDS  LENGTH  OF  SIDE  OF  CUBE  WITH  CENTER  <H,K,L)  WHICH  CONTAINS  P3  OF 
11740  i  THE  NORMAL  DISTRIBUTION  WITH  MEAN  0  AND  STANDARD  DEVIATIONS  <S1,S2,S3). 
11745  !  FOR  FNEllrc  CHECK  T6.  IF  T6>40  RESULT  SUSPECT. 

11750  !  SUBROUTINES  NEEDED:  ERF,  ERF3,  ERFC1. 

11755  DEF  FNE1 1 rc <H, K , L , S 1 , S2 , S3 , P3 , Run  n, Rmax , T6 > 

11760  T6=0 

11765  IF  < P3  >  1  )  OR  < P 3 < 0 )  THEN  RETURN  -1E99 

11770  IF  P3= 1  THEN  RETURN  1E99 

11775  IF  P3=0  THEN  RETURN  0 

11780  Sq=l . 41421356237 

11785  Ul-10Sq*Sl) 
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11790  U2=l/(Sq*S2> 

11795  Sqpi =2. 50662827463 
11800  U3=l/<Sq*S3> 

11805  Y 1 = 1 / ( Sqp i  *S  1  > 

11810  Y2= 1 / <  Sqp i  *S2 ) 

11815  Y3= 1 / <  Sqp i *S3 ) 

11820  Rl=Rmax 
11825  R0=Rmin 
11830  T6=0 

11835  A=.5*<R1+R0) 

11840  GOSUB  11930 

11845  IF  CABS  (UX.  1  *P3  >  OR  <  BBS  <  R  1  -R0  >  <  .  0 1  *R0  >  THEN  11865 
11850  GOTO  11835 
11855  A=. 5*<R1+R0 > 

11860  GOSUB  11930 

11865  IF  C0  =  0  THEN  C0  =  EXP < -2*H*Ax < S 1 *S  1  > > 

11870  IF  C2=0  THEN  C2=EXP < -2*K*A/ < S2*S2 > ) 

11875  IF  C4=0  THEN  C4=EXP < -2*L*A/ < S3*S3 ) > 

11880  B1=Y1*C1*< 1+C0)*F2*F3+Y2*C3*< 1 +C2 > *F 1 *F3  + Y3*C5* ( 1 +C4 > *F 1 *F2 

11885  IF  B 1 <  =0  THEN  11835 
11890  A1=A-U/B1 

11895  IF  ABS<F-P3>< IE-10  THEN  RETURN  A1 

11900  IF  <ABS(A1-AX5E-5*A1>  AND  <ABS(F-P3X 1E-3*P3>  THEN  RETURN  A1 
11905  IF  ABS<Al-AX  =  lE-8  THEN  RETURN  A1 
11910  A=A 1 

11915  IF  T6  >40  THEN  RETURN  A1 

11920  IF  ( A<  R0 )  OR  <  A  >R 1 >  THEN  11855 

11925  GOTO  11860 

11930  T0= ( H+A ) *U 1 

11935  T 1 = ( H- A ) *U 1 

11940  T2= ( K+A  >  *U2 

11945  T3=  <  K-A ) *U2 

11950  T4= ( L+ A ) *U3 

11955  T5= ( L- A ) *U3 

11960  C1=F  =  EXP<-T1*T1  > 

11965  C3=F=EXP<-T3*T3> 

11970  C5=F=EXP(-T5*T5> 

11975  C0=C2=C4=0 

11980  IF  T 1 < 0  THEN  12005 
11985  IF  C1=0  THEN  GOSUB  12085 
1  1990  C0  =  EXP<-2*H*A/<S1*S1  ) ) 

1  1995  Fl  =  . 5*Cl*<FNErfc 1 (T1 , 1 >-C0*FNErfc 1 <T0,  1  >  > 

12000  GOTO  12010 

12005  Fl=. 5*<FNErf <T0)-FNErf 3<T1 , Cl  >  > 

12010  IF  T3<  0  THEN  12035 

12015  IF  C3=0  THEN  GOSUB  12085 

12020  C2  =  EXP<-2*K*AXS2*S2)  > 

12025  F2=.5*C3*(FNErfcl(T3, 1  >-C2*FNErfc 1 <T2, 1 >  > 

12030  GOTO  12040 

12035  F2=.  5*<FNErf <T2>-FNErf 3<T3, C3>  > 

12040  IF  T5<  0  THEN  12065 

12045  IF  C5=0  THEN  GOSUB  12085 

12050  C4  =  EXP(-2*L*AXS3*S3>  ) 

12055  F3*. 5*C5*<FNEr fc 1 CT5, 1 >-C4*FNErfc 1 <T4, 1 >  > 
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12060  GOTO  12070 

12065  F3=.5*(FNErf(T4>-FNErf3(T5,C5>) 

12070  F=F1 *F2*F3 

12075  U=F-P3 

12080  IF  U  >0  THEN  12095 

12085  R0= A 

12090  GOTO  1210O 

12095  R 1 =A 

12100  T6=T6+ 1 

12105  IF  ABSCF-. 5>=. 5  THEN  11835 
12110  RETURN 
12115  FNEND 

12120  !  ***********FNXerfc<X>  =  X*EXP(X*X>*ERFC<X)************ 

12125  DEF  FNXerf  c  <  X  > 

12130  X2=X*X 

12135  IF  ABS ( X  )  > .  5  THEN  12150 

12140  Erfc 1  =  . 5-X*<  <  03. 560984370 18E-2*X2 rg . 99638348862 > *X2  +  2 1 . 9792616183>*X2+2 

42.667955231  >/< << X2+ 1 5 . 0827976304 > *X2+9 1 . 1 64  9054045 > *X2  +  2 1 5 . 05887587 >  +  . 5 
12145  RETURN  X*EXP < X2 ) *Er f c  1 
12150  K4=ABS(X) 

12155  IF  K4  >4  THEN  12180 

12160  K3=<  < <  < -6. 0858 151 9597E-6*K4+. 56437 1 606864 > *K4+4 . 2677201 0709 >*K4+ 14. 57189 

85969>*K4+26. 0947469561 )  *  K  4  +  2  2 . 89^9928517 

12165  Erfc 1=K3'C  < <  <  < K4  +  7 . 56884822936 ) *K4  +  26 . 2887957588 > *K4  +  50 . 2732028638 > *K4  +  5 

1. 9335706876 >*K4+22. 8989857499) 

12170  IF  X<0  THEN  RETURN  X* < 2*EXP ( X2 ) -Erf c 1  ) 

12175  RETURN  X*Erf c  1 

12180  IF  X< -4 . 900  THEN  RETURN  2*X*EXP<X*X> 

12185  K6=l/X2 

12190  Bl  =  . 564189583548  !  1/SGR<PI) 

12195  Erfc 1=<B1-K6*< ( (3. 243 1 95 1 9278E-2*K6+ . 2439 1 1 029489 ) *K6+ . 1 1 9903955268 > *K6+ 

•  012130827639)/<  <  (  K6+  1 . 4377  1  227937  )  *K6  +  .  48955244  1  96  1  >*K6  +  4. 3002664  3453E-2>>--'K4 
12200  IF  X<0  THEN  RETURN  X* < 2*EXP < X2 > -Er f c  1  > 

12205  RETURN  X*Erfc 1 
12210  FNEND 

12215  !  ********FNDxerf3<X, Exp>  =  ERF<X>/X;  EXPOX*X>  AVAILABLE  IF  Exp  >  0.** 

12220  DEF  FNDxerf 3<X, Exp) 

12225  X2=X*X 

12230  IF  ABS  <  X  >  > . 5  THEN  12245 

12235  Erf=<  <  < -3 . 560984370 1 8E-2*X2  +  6 . 99638348862 > *X2  +  2 1 .9792616i83>*X2+242.6679 

55231 > /<<(X2+15. 0827976304 )*X2+91. 1 649054045 > *X2+2 1 5 . 05887587 > 

12240  RETURN  Erf 

12245  K4  =  ABS  <  X ) 

12250  IF  K4  >4  THEN  12275 

12255  K3=  < ( < (-6. 0858 1 5 1 9597E-6*K4+ . 564371 6068 64 >*K4  +  4. 267720 10709 >*K4  + 14. 57189 

85969 )*K4+26. 0947469561 )*K4+22. 8989928517 

12260  K3  =  K3/<  <  < ( < K4  +  7 . 5688  4822936 > *K4  +  26 . 2887957588 > *K4  +  50 . 2732028638 > *K 4  +  5 1 . 9 

335706876 )*K4+22. 8989857499) 

12265  IF  Exp<  0  THEN  Exp  =  EXPOX2> 

12270  RETURN  < 1 -Exp*K3 > xK4 
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12275  IF  K4  >  =  4 . 883  THEN  RETURN  1/K4 
12280  K6*l/X2 

12285  Bl*. 564189583548 !  1/SQR(PI) 

12290  Erf =  (B1 -K6*(  (  <3. 243 1951 9278E-2*K6+. 24391 1029489 )*K6+. 1 1 9903955268 ) *K6  + . 0 

121 30827639 ) /(  (  (K6+1 . 4377 1 227937 ) «K6  + . 48955244 1 96 1 ) *K6  +  4 . 30026643453E-2 ) )/K4 
12295  IF  Exp<  0  THEN  Exp=EXP(-X2) 

12300  RETURN  ( 1 -Exp*Erf ) /K4 

12305  FNEND 

12310  !  **********#****DXD8WS  <X)  *****■*•*■*■******♦■*■ 

12315  !  Daws ( X )  =  1 /X*EXP ( -X*X ) *  I NTEGRRL  FROM  0  TO  X  OF  EXP(T*T)  DT. 

12316  !  REF:  MATH  OF  COMP.  1970,  PP.  171-178.  DAWSON ' S  INTEGRAL/X 
12320  DEF  FNDxdaws(Kl) 

12325  Daws=Kl*Kl 

12333  IF  ABS C K 1 ) >2 . 5  THEN  12350 

12335  K4“ ( ( << -2. 38594565696E-2*Daws+l .677951 1 6 1 89 ) »Daws- 1 5 . 1 982 1 52422 ) *Daws+4 1 

9. 67290228 )*Daws- 1284. 05832279 )*Daws+ 10832. 6558873 

12340  Daws  =  K4/<  < <  <  < Daws+ 1 9 . 9422336364 ) *Daws  +  2 1 9 . 72833 1 833 > *Daws+ 1 489 . 43557242 > 

♦Daws+5937. 7 1276935 >*Daws+ 10832. 6558 772 > 

12345  RETURN  Daws 

12350  IF  ABS < K 1 >  >3 . 5  THEN  12380 

12355  K4=88. 76 1 9386764/ <Daws+3. 47393742586) 

12360  K4= 13. 8216341 182/ (Daws- 12. 68 17901 598+K4) 

12365  K4=37. 39027635 1 2/ (Daws+4. 07068 101667+K4) 

12370  Daws=( . 2 1 422 1 965778/ < Daws-4 . 9 1 60536574 1 +K4 )+ . 500652754437 >/Daus 

12375  RETURN  Daws 

12380  IF  ABS ( K 1 ) >=5  THEN  12405 

12385  K4=17. 76 1 778077/ ( Daws-4 . 1083233 770) 

12390  K4=-. 53714272981 /(Daws- 10. 385024821 +K4) 

12395  Daws=( . 2492463242 1 /( Daws- 1 . 5853935006+K4 )+. 50000965 1 99 ) /Daws 

12400  RETURN  Daws 

12405  K4=-4. 1 2544065608/ ( Daws- 1 1 . 1952164237) 

12410  K4=-2. 48787658804/ (Daws-4. 6731 2022141 +K4) 

12415  Daws  =  .  5*(  l  +  (  .  749999 1 90567/ ( Daws-2 . 500 1 7 1 1 6686  +  K4 )  +  . 50000000 1 674 ) /Daws ) /D 

aws 

12420  RETURN  Daws 
12425  FNEND 
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